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(nullptr), _firstEntryOfCurrentFile(-1), _detInfo(nullptr) {
  _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 = nullptr;
  _detInfo = nullptr;
  _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 (_CheckBranch(detBranchName)) {

    int detIndex = GetDetectorIndex(detector); // detIndex == -1 if the detector has not been set
    TBranch *partHitsBranch = _chain->GetBranch(std::string(detector).append("-GGSPartHits").c_str());
    TBranch *posHitsBranch = _chain->GetBranch(std::string(detector).append("-GGSPosHits").c_str());

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

    // Position hits
    if (readPosHits) {
      if (posHitsBranch) {
        infoString.append(" (+PosHits)");
        if (detIndex == -1) {
          _posHitArrays.push_back(nullptr);
          posHitsBranch->SetAddress(&(_posHitArrays.back()));
          _readPosHits.push_back(true);
        }
        else {
          posHitsBranch->SetAddress(&(_posHitArrays[detIndex]));
          _readPosHits[detIndex] = true;
        }
        branches[2] = posHitsBranch;
      }
      else {
        throw std::runtime_error(std::string("PosHits missing for detector ") + detector);
      }
    }
    else {
      if (detIndex == -1) {
        _posHitArrays.push_back(nullptr);
        _readPosHits.push_back(false);
      }
      else {
        if (posHitsBranch) {
          posHitsBranch->ResetAddress();
        }
        _posHitArrays[detIndex] = nullptr; // FIXME: should the element be deleted before this? Is it owned by the file/branch?
        _readPosHits[detIndex] = false;
      }
    }

    // Particle hits
    if (readPartHits) {
      if (partHitsBranch) {
        infoString.append(" (+PartHits)");
        if (detIndex == -1) {
          _partHitArrays.push_back(nullptr);
          partHitsBranch->SetAddress(&(_partHitArrays.back()));
          _readPartHits.push_back(true);
        }
        else {
          partHitsBranch->SetAddress(&(_partHitArrays[detIndex]));
          _readPartHits[detIndex] = true;
        }
        branches[1] = partHitsBranch;
      }
      else {
        throw std::runtime_error(std::string("PartHits missing for detector ") + detector);
      }
    }
    else {
      if (detIndex == -1) {
        _partHitArrays.push_back(nullptr);
        _readPartHits.push_back(false);
      }
      else {
        if (partHitsBranch) {
          partHitsBranch->ResetAddress();
        }
        _partHitArrays[detIndex] = nullptr; // FIXME: should the element be deleted before this? Is it owned by the file/branch?
        _readPartHits[detIndex] = false;
      }
    }

    if (detIndex == -1) {
      _intHitArrays.push_back(nullptr);
      _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));
      GGSSmartLog::MuteOutput();
      bool isDetBranch = _CheckBranch(branch->GetName());
      GGSSmartLog::UnmuteOutput();
      if (isDetBranch) {
        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 nullptr;
}

GGSTIntHit* GGSTHitsReader::GetHit(int iDet, unsigned int iHit) {
  static TClonesArray *hitArray = nullptr;
  hitArray = _GetHitArray(iDet);
  if (hitArray) {
    static GGSTIntHit *hit = nullptr;
    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 = nullptr;
      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 nullptr;
}

//....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 = nullptr;
      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 = nullptr;
  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.....

TClonesArray* GGSTHitsReader::_GetHitArray(int iHitArray) {

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

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

GGSTHitDetInfo* GGSTHitsReader::_GetDetectorInfo(int iDet) {

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