GGSTHitsReader.cpp

Go to the documentation of this file.

/*
 * GGSTHitsReader.cpp
 *
 *  Created on: Nov 8, 2010
 *      Author: Elena Vannuccini
 *  Reworked on: 17 Aug 2011
 *      Author: Nicola Mori
 *      Content: MC truth stripped out. File opening removed.
 */

#include "utils/GGSSmartLog.h"
#include "montecarlo/readers/GGSTHitsReader.h"

#include "TFile.h"
#include "TFriendElement.h"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

GGSTHitsReader::GGSTHitsReader() :
    _chain(NULL), _firstEntryOfCurrentFile(-1), _detInfo(NULL) {
  _detectors.reserve(50);
  // The elements of these vectors will be used as output buffers for branches, so
  // their location in memory must not change. The instructions below guarantee
  // correct behavior for a number of detectors up to 50.
  _intHitArrays.reserve(50);
  _partHitArrays.reserve(50);
  _posHitArrays.reserve(50);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

GGSTHitsReader::~GGSTHitsReader() {

  for (unsigned int i = 0; i < _intHitArrays.size(); i++)
    delete _intHitArrays[i];
  for (unsigned int i = 0; i < _partHitArrays.size(); i++)
    delete _partHitArrays[i];
  for (unsigned int i = 0; i < _posHitArrays.size(); i++)
    delete _posHitArrays[i];
  delete _detInfo;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

bool GGSTHitsReader::SetChain(TChain *hitsChain) {

  _chain = hitsChain;
  _detInfo = (TClonesArray*) (_chain->GetFile()->Get("GGSHitDetInfo"));

  // Search for hit branches in master chain
  for (int iBranch = 0; iBranch < _chain->GetListOfBranches()->GetEntries(); iBranch++) {
    if (TString(((TBranch*) (_chain->GetListOfBranches()->At(iBranch)))->GetName()).Contains("GGSIntHits")) {
      return true;
    }
  }

  // Search for hit branches in friend chains
  if (_chain->GetListOfFriends()) {
    for (int iFriend = 0; iFriend < _chain->GetListOfFriends()->GetEntries(); iFriend++) {
      TTree *currFriend = ((TFriendElement*) (_chain->GetListOfFriends()->At(iFriend)))->GetTree();
      for (int iBranch = 0; iBranch < currFriend->GetListOfBranches()->GetEntries(); iBranch++) {
        if (TString(((TBranch*) (currFriend->GetListOfBranches()->At(iBranch)))->GetName()).Contains("GGSIntHits")) {
          return true;
        }
      }
    }
  }

  _chain = NULL;
  _detInfo = NULL;
  _firstEntryOfCurrentFile = -1;

  return false;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

void GGSTHitsReader::GetEntry(Long64_t entry) {
  static const std::string routineName("GGSTHitsReader::GetEntry");

  if (entry != _chain->GetReadEntry()) {
    _chain->GetEntry(entry);
  }

  if (_chain->GetChainEntryNumber(0) != _firstEntryOfCurrentFile) {
    COUT(DEBUG) << "Switching to file " << _chain->GetFile()->GetName() << ENDL;
    _firstEntryOfCurrentFile = _chain->GetChainEntryNumber(0);

    // Read the GGSHitDetInfo TClonesArray from the new file
    _detInfo = (TClonesArray*) (_chain->GetFile()->Get("GGSHitDetInfo"));

    // Update the detector indexes in _detectors
    for (unsigned int iDet = 0; iDet < _detectors.size(); iDet++)
      _detectors[iDet].detInfoArrayIndex = _GetDetectorIndexFromStorage(_detectors[iDet].logVolName.data());

    // Reload time bins and thresholds
    _timeBins.clear();
    _thresholds.clear();
    for (unsigned int iDet = 0; iDet < _detectors.size(); iDet++) {
      delete _timeBins[iDet];
      _timeBins.push_back((TArrayF*) (_chain->GetFile()->Get(TString(_detectors[iDet].name) + "-GGSIntHits-TimeBins")));
      delete _thresholds[iDet];
      _thresholds.push_back(
          (std::vector<float>*) (_chain->GetFile()->Get(TString(_detectors[iDet].name) + "-GGSIntHits-Thresholds")));
    }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

void GGSTHitsReader::SetDetector(const char* detector, bool readPartHits, bool readPosHits) {

  static const std::string routineName("GGSTHitsReader::SetDetector");

  if (!_chain)
    return;

  if (readPosHits) {
    readPartHits = true;
  }

  TString detBranchName = detector;
  detBranchName.Append("-GGSIntHits");
  TBranch *branches[3] = { nullptr, nullptr, nullptr }; // New branches (for integrated, particle and position hits)
  if (_CheckDetector(detector) && _CheckBranch(detBranchName)) {

    std::string infoString("Set detector ");
    infoString.append(detector);

    if (readPartHits) {
      TString branchName(detector);
      branchName.Append("-GGSPartHits");
      TBranch *branch = _chain->GetBranch(branchName);
      if (branch) {
        infoString.append(" (+PartHits)");
        _partHitArrays.push_back(NULL);
        _chain->SetBranchAddress(branchName, &(_partHitArrays.back()));
        branches[1] = _chain->GetBranch(branchName);
        _readPartHits.push_back(true);
      }
      else {
        throw std::runtime_error(std::string("PartHits missing for detector ") + detector);
      }
    }
    else {
      _partHitArrays.push_back(NULL);
      _readPartHits.push_back(false);
    }

    if (readPosHits) {
      TString branchName(detector);
      branchName.Append("-GGSPosHits");
      TBranch *branch = _chain->GetBranch(branchName);
      if (branch) {
        infoString.append(" (+PosHits)");
        _posHitArrays.push_back(NULL);
        _chain->SetBranchAddress(branchName, &(_posHitArrays.back()));
        branches[2] = _chain->GetBranch(branchName);
        _readPosHits.push_back(true);
      }
      else {
        _partHitArrays.pop_back();
        _readPartHits.pop_back();
        _chain->ResetBranchAddress(_chain->GetBranch((std::string(detector) + "-GGSPartHits").c_str()));
        throw std::runtime_error(std::string("PosHits missing for detector ") + detector);
      }
    }
    else {
      _posHitArrays.push_back(NULL);
      _readPosHits.push_back(false);
    }

    _intHitArrays.push_back(NULL);
    _chain->SetBranchAddress(detBranchName, &(_intHitArrays.back()));
    branches[0] = _chain->GetBranch(detBranchName);
    _detInfo = (TClonesArray*) (_chain->GetFile()->Get("GGSHitDetInfo"));
    _timeBins.push_back((TArrayF*) (_chain->GetFile()->Get(TString(detector) + "-GGSIntHits-TimeBins")));
    _thresholds.push_back((std::vector<float>*) (_chain->GetFile()->Get(TString(detector) + "-GGSIntHits-Thresholds")));

    Detector det;
    det.name = detector;
    det.logVolName = detector;
    det.logVolName = det.logVolName.substr(0, det.logVolName.find_first_of('.')); // remove names of scorer and sensitive detector
    det.detInfoArrayIndex = _GetDetectorIndexFromStorage(det.logVolName.c_str());
    _detectors.push_back(det);

    COUT(INFO) << infoString << ENDL;
  }
  else {
    throw std::runtime_error(std::string("IntHits missing for detector ") + detector);
  }

  // Synchronize the branches that have just been set for reading
  if (_chain->GetReadEntry() != -1) {
    for (int iBranch = 0; iBranch < 3; ++iBranch) {
      if (branches[iBranch]) {
        branches[iBranch]->GetEntry(_chain->GetReadEntry());
      }
    }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

bool GGSTHitsReader::DetectorExists(const char *detector) {

  if (_chain) {
    TString detectorString(detector);
    detectorString.Append("-GGSIntHits");
    return _CheckBranch(detectorString.Data());
  }

  return false;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

void GGSTHitsReader::ListDetectors() {

  static const std::string routineName("GGSTHitsReader::ListDetectors");

  if (_chain) {
    TObjArray *branches = _chain->GetListOfBranches();
    COUT(INFO) << "List of detectors in hits tree:" << ENDL;
    for (int i = 0; i < branches->GetEntries(); i++) {
      TBranch *branch = (TBranch*) (branches->At(i));
      if (_CheckBranch(branch->GetName())) {
        TString detName = branch->GetName();
        detName.Remove(detName.Length() - 11); // Remove the trailing "-GGSIntHits"
        std::string infoString("   ");
        infoString.append(detName.Data());
        if (_chain->GetBranch(detName + "-GGSPartHits"))
          infoString.append("  (PartHits)");
        if (_chain->GetBranch(detName + "-GGSPosHits"))
          infoString.append("  (PosHits)");
        CCOUT(INFO) << infoString << ENDL;
      }
    }

  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

std::vector<TString> GGSTHitsReader::GetListOfDetectors() {

  static const std::string routineName("GGSTHitsReader::GetListOfDetectors");

  std::vector<TString> _list;

  if (_chain) {
    TObjArray *branches = _chain->GetListOfBranches();
    for (int i = 0; i < branches->GetEntries(); i++) {
      TBranch *branch = (TBranch*) (branches->At(i));
      if (_CheckBranch(branch->GetName())) {
        TString detName = branch->GetName();
        detName.Remove(detName.Length() - 11); // Remove the trailing "-GGSIntHits"
        _list.emplace_back(detName);
      }
    }

  }

  return _list;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

Int_t GGSTHitsReader::GetNHits(const char *detector) {
  static short int iHitArray;
  iHitArray = GetDetectorIndex(detector);
  if (iHitArray != -1)
    return (_GetHitArray(iHitArray))->GetEntries();
  else
    return 0;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

Int_t GGSTHitsReader::GetNHits(int iDet) {

  TClonesArray *detArray = _GetHitArray(iDet);
  if (detArray)
    return detArray->GetEntries();
  else
    return 0;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

GGSTIntHit *GGSTHitsReader::GetHit(const char *detector, unsigned int iHit) {
  int iDet = GetDetectorIndex(detector);
  if (iDet != -1)
    return GetHit(iDet, iHit);
  else
    return NULL;
}

GGSTIntHit *GGSTHitsReader::GetHit(int iDet, unsigned int iHit) {
  static TClonesArray* hitArray = NULL;
  hitArray = _GetHitArray(iDet);
  if (hitArray) {
    static GGSTIntHit *hit = NULL;
    hit = (GGSTIntHit*) (hitArray->At(iHit));
    if (hit) {
      // Compute non-persistent informations
      hit->eDep = hit->eDepTimeBin.GetSum();
      hit->_hitDetInfo = _GetDetectorInfo(_detectors[iDet].detInfoArrayIndex);
      hit->_timeBins = _timeBins[iDet];
      if (hit->_hitDetInfo)
        hit->_hitVolInfo = (GGSTHitVolInfo*) (hit->_hitDetInfo->volumes.At(hit->_volumeIndex));
      else
        hit->_hitVolInfo = NULL;
      if (_readPartHits[iDet]) {
        hit->_partHits = (TClonesArray*) (_partHitArrays[iDet]->At(iHit));
      }
      if (_readPosHits[iDet]) {
        // Set pointer to position hits collection for each particle hit
        TClonesArray *partHits = hit->_partHits;
        int nPartHits = partHits->GetEntries();
        for (int iPartHit = 0; iPartHit < nPartHits; iPartHit++) {
          GGSTPartHit *partHit = (GGSTPartHit*) (partHits->At(iPartHit));
          partHit->_posHits = (TClonesArray*) (_posHitArrays[iDet]->At(partHit->_posHitIndex));
        }
      }
      return hit;
    }
  }
  return NULL;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

int GGSTHitsReader::GetDetectorIndex(const char* detector) {
  for (UInt_t i = 0; i < _detectors.size(); i++) {
    if (!strcmp(_detectors[i].name.data(), detector))
      return i;
  }
  return -1;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

bool GGSTHitsReader::HasPartHits(const char* detector) {
  if (_chain) {
    TString branchName(detector);
    branchName.Append("-GGSPartHits");
    TBranch *branch = _chain->GetBranch(branchName);
    if (branch)
      return true;
    else
      return false;
  }
  else
    return false;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

bool GGSTHitsReader::HasPosHits(const char* detector) {
  if (_chain) {
    TString branchName(detector);
    branchName.Append("-GGSPosHits");
    TBranch *branch = _chain->GetBranch(branchName);
    if (branch)
      return true;
    else
      return false;
  }
  else
    return false;
}

float GGSTHitsReader::GetIntHitThreshold(const std::string &detector) {
  const std::string routineName("GSTHitsReader::GetIntHitThreshold");
  int detIndex = GetDetectorIndex(detector.c_str());
  if (detIndex == -1) {
    throw std::range_error(routineName + ": detector " + detector + " not found");
  }

  if (_thresholds[detIndex]) {
    return (*(_thresholds[detIndex]))[0];
  }
  else {
    throw std::runtime_error(routineName + ": thresholds for integrated hits not available for detector " + detector);
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

float GGSTHitsReader::GetPartHitThreshold(const std::string &detector) {
  const std::string routineName("GSTHitsReader::GetIntHitThreshold");
  int detIndex = GetDetectorIndex(detector.c_str());
  if (detIndex == -1) {
    throw std::range_error(routineName + ": detector " + detector + " not found");
  }

  if (_thresholds[detIndex]) {
    return (*(_thresholds[detIndex]))[1];
  }
  else {
    throw std::runtime_error(routineName + ": thresholds for particle hits not available for detector " + detector);
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

float GGSTHitsReader::GetPosHitThreshold(const std::string &detector) {
  const std::string routineName("GSTHitsReader::GetIntHitThreshold");
  int detIndex = GetDetectorIndex(detector.c_str());
  if (detIndex == -1) {
    throw std::range_error(routineName + ": detector " + detector + " not found");
  }

  if (_thresholds[detIndex]) {
    return (*(_thresholds[detIndex]))[2];
  }
  else {
    throw std::runtime_error(routineName + ": thresholds for position hits not available for detector " + detector);
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

int GGSTHitsReader::_GetDetectorIndexFromStorage(const char* detector) {

  if (_detInfo) {
    for (int iDet = 0; iDet < _detInfo->GetEntries(); iDet++) {
      static GGSTHitDetInfo *det = NULL;
      det = (GGSTHitDetInfo*) (_detInfo->At(iDet));
      if (det && det->detectorName == detector)
        return iDet;
    }
  }
  return -1;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

bool GGSTHitsReader::_CheckBranch(const char *branchName) {

  static const std::string routineName("GGSTHitsReader::_CheckBranch");
  TBranch* branch = _chain->GetBranch(branchName);

  // Check whether the branch exists or not
  if (!branch) {
    //COUT << routineName << "Warning: branch " << branchName << " not found." << ENDL;
    return false;
  }

  // Check whether the branch is a TClonesArray or not
  if (strcmp(branch->GetClassName(), "TClonesArray")) {
    //COUT << routineName << "Warning: branch " << branchName << " is not a valid hit container (TClonesArray)." << ENDL;
    return false;
  }

  // Check whether the TClonesArray contains hit objects or not
  char *oldBuffer = branch->GetAddress();
  Long64_t oldReadEntry = branch->GetReadEntry();
  TClonesArray *buffer = NULL;
  branch->SetAddress(&buffer);
  branch->GetEntry(0);
  branch->SetAddress(oldBuffer);
  branch->GetEntry(oldReadEntry);
  std::string hitClassName(buffer->GetName());
  hitClassName.erase(hitClassName.size() - 1); // Remove the trailing 's' that TClonesArray name places after the class name
  TClass *clonesClass = TClass::GetClass(hitClassName.c_str());
  if (!clonesClass) {
    COUT(ERROR) << "No dictionary for class " << hitClassName << " inside the TCloneArray in branch " << branchName
        << ENDL;
    return false;
  }
  if (!(clonesClass->GetBaseClass("GGSTIntHit"))) {
    COUT(ERROR) << "Class " << hitClassName << " inside the TCloneArray in branch " << branchName
        << " does not inherit from GGSTIntHit" << ENDL;
    return false;
  }
  return true;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

bool GGSTHitsReader::_CheckDetector(const char *detector) {

// Check whether the detector has already been set or not
  for (unsigned int i = 0; i < _detectors.size(); i++) {
    if (_detectors[i].name == detector)
      return false; // The branch has not to be added again.
  }

  return true;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

TClonesArray* GGSTHitsReader::_GetHitArray(int iHitArray) {

  if (iHitArray > (int) (_intHitArrays.size() - 1) || iHitArray < 0)
    return NULL;
  return _intHitArrays[iHitArray];
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

GGSTHitDetInfo *GGSTHitsReader::_GetDetectorInfo(int iDet) {

  if (_detInfo) {
    static GGSTHitDetInfo *det = NULL;
    det = (GGSTHitDetInfo*) (_detInfo->At(iDet));
    return det;
  }
  return NULL;
}