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 "montecarlo/readers/GGSTHitsReader.h"

#include "TFile.h"
#include "TFriendElement.h"
#include "utils/GGSSmartLog.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;
  delete _detInfo;
  _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;
  delete _detInfo;
  _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) {
    GGSCOUT(DEBUG) << "Switching to file " << _chain->GetFile()->GetName() << GGSENDL;
    _firstEntryOfCurrentFile = _chain->GetChainEntryNumber(0);

    // Read the GGSHitDetInfo TClonesArray from the new file
    delete _detInfo;
    _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) == 0) {
    int detIndex = GetDetectorIndex(detector); // detIndex == -1 if the detector has not been set
    std::string partHitsBranchName{std::string(detector).append("-GGSPartHits")};
    std::string posHitsBranchName{std::string(detector).append("-GGSPartHits")};
    TBranch *partHitsBranch = _chain->GetBranch(partHitsBranchName.c_str());
    TBranch *posHitsBranch = _chain->GetBranch(posHitsBranchName.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);
          _chain->SetBranchAddress(posHitsBranchName.c_str(), &(_posHitArrays.back()));
          _readPosHits.push_back(true);
        } else {
          _chain->SetBranchAddress(posHitsBranchName.c_str(), &(_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);
          _chain->SetBranchAddress(partHitsBranchName.c_str(), &(_partHitArrays.back()));
          _readPartHits.push_back(true);
        } else {
          _chain->SetBranchAddress(partHitsBranchName.c_str(), &(_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);
      delete _detInfo;
      _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);
    }
    GGSCOUT(INFO) << infoString << GGSENDL;
  } 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, bool particleHits, bool positionHits) {
  int type = -1;
  if (_chain) {
    std::string detectorString{detector};
    detectorString += "-GGSIntHits";
    type = _CheckBranch(detectorString.c_str());
  }
  if (type == 0 && particleHits) {
    std::string detectorString{detector};
    detectorString += "-GGSPartHits";
    type = _CheckBranch(detectorString.c_str());
  }
  if (type == 1 && positionHits) {
    std::string detectorString{detector};
    detectorString += "-GGSPosHits";
    type = _CheckBranch(detectorString.c_str());
  }

  return (type != -1);
}

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

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

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

//....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));
      std::string branchName(branch->GetName());
      if (_CheckBranch(branch->GetName()) == 0) {
        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.....

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

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

std::string GGSTHitsReader::GetIntHitsClassName(const char *detector) {
  int iDet = GetDetectorIndex(detector);
  if (iDet >= 0) {
    return GetIntHitsClassName(iDet);
  } else {
    return "";
  }
}

std::string GGSTHitsReader::GetIntHitsClassName(int iDet) {
  TClonesArray *hitArray = _GetHitArray(iDet);
  if (hitArray) {
    std::string retValue = hitArray->GetName();
    return retValue.substr(0, retValue.size() - 1);
  } else {
    return "";
  }
}
//....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.....

int GGSTHitsReader::_CheckBranch(const char *branchName) {
  const std::string routineName("GGSTHitsReader::_CheckBranch");

  std::string branchNameStr{branchName};
  int type = -1; // 0: integrated, 1: particle, 2: position
  if (branchNameStr.find("-GGSIntHits") != std::string::npos) {
    type = 0;

  } else if (branchNameStr.find("-GGSPartHits") != std::string::npos) {
    type = 1;
  } else if (branchNameStr.find("-GGSPosHits") != std::string::npos) {
    type = 2;
  } else {
    return -1;
  }

  TBranch *branch = _chain->GetBranch(branchName);

  // Check whether the branch exists or not
  if (!branch) {
    // GGSCOUT << routineName << "Warning: branch " << branchName << " not found."
    // << GGSENDL;
    return -1;
  }

  // Check whether the branch is a TClonesArray or not
  if (type == 0) {
    if (strcmp(branch->GetClassName(), "TClonesArray")) {
      return -1;
    }
  } else {
    if (strcmp(branch->GetClassName(), "TObjArray")) {
      return -1;
    }
  }

  if (type == 0) {
    // Check whether the TClonesArray contains hit objects or not
    TClonesArray **bufferAddr =
        (TClonesArray **)(branch->GetAddress()); // See the documentation of TBranchElement::GetAddress
    bool addressNotSet = false;
    if (!bufferAddr) {
      bufferAddr = new TClonesArray *;
      *bufferAddr = nullptr;
      branch->SetAddress(bufferAddr);
      addressNotSet = true;
    }

    auto oldReadEntry = branch->GetReadEntry();
    branch->GetEntry(0);
    std::string hitClassName((*bufferAddr)->GetName());
    if (addressNotSet) {
      branch->ResetAddress();
      bufferAddr = nullptr;
    } else {
      branch->GetEntry(oldReadEntry);
      assert(branch->GetReadEntry() == oldReadEntry);
    }
    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) {
      GGSCOUT(ERROR) << "No dictionary for class " << hitClassName << " inside the TCloneArray in branch " << branchName
                     << GGSENDL;
      return -1;
    }
    if (!(clonesClass->GetBaseClass("GGSTIntHit"))) {
      GGSCOUT(ERROR) << "Class " << hitClassName << " inside the TCloneArray in branch " << branchName
                     << " does not inherit from GGSTIntHit" << GGSENDL;
      return -1;
    }
  }

  return type;
}

//....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;
}