Dodawanie własnego detektora

Dodawanie własnego detektora część I

Pierwszym krokiem w konstrukcji własnego detektora jest stworzenie własnej biblioteki (patrz poprzednie lekcje). Tutaj tworzymy klasy dziedziczące po:

FairMCPoint do przechowywania "punktów MC" z naszego detektora
FairParGenericSet do przechowywania parametrów (w zasadzie nie jest ona używana ale FairROOT robi bez niej problemy)
FairGeoSet klasa reprezentująca geometrię detektora
FairContFact klasa produkująca kontenery dla naszego detektora
FairDetector właściwa klasa reprezentująca detektor

Pierwsze 4 klasy w zasadzie są wymagane przez framework i można je skopiować z innego detektora "inteligetnie" zmieniając nazwy parametrów itd. Z naszego punktu widzenia tak naprawdę interesuje nas to co dziedziczy po FairDetector. Tu na tapetę weźmiemy klasę MpdMcord z eksperymentu MPD. Header tej klasy wygląda tak:

#ifndef MCORD_MCORD_MPDMCORD_H_
#define MCORD_MCORD_MPDMCORD_H_
#include "FairDetector.h"
#include "TClonesArray.h"
#include "TLorentzVector.h"
#include "TVector3.h"

class MpdMCordHit;
class FairVolume;
class MpdMcordGeo;
class MpdMcordPoint;

class MpdMcord : public FairDetector {
  Int_t fTrackID;       //!  track index
  Int_t fVolumeID;      //!  volume id
  TLorentzVector fPos;  //!  position
  TLorentzVector fMom;  //!  momentum
  Double_t fTime;       //!  time
  Double_t fLength;     //!  length
  Double_t fELoss;      //!  energy loss

public:
  MpdMcord();
  MpdMcord(const char* name, Bool_t active);
  virtual ~MpdMcord();
  // Defines the action to be taken when a step is inside the
  // active volume. Creates MpdStsPoints and adds them to the collection.
  // @param vol  Pointer to the active volume
  virtual Bool_t ProcessHits(FairVolume* vol = 0);

  // If verbosity level is set, print hit collection at the
  // end of the event and resets it afterwards.
  virtual void EndOfEvent();

  // Registers the hit collection in the ROOT manager.
  virtual void Register();

  // Accessor to the hit collection
  virtual TClonesArray* GetCollection(Int_t iColl) const;

  // Screen output of hit collection.
  virtual void Print() const;

  // Clears the hit collection
  virtual void Reset();

  // *@param cl1     Origin
  // *@param cl2     Target
  // *@param offset  Index offset
  virtual void CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset);

  // Constructs the STS geometry
  virtual void ConstructGeometry();


  virtual Bool_t CheckIfSensitive(std::string name);
  MpdMcordPoint* AddHit(Int_t trackID, Int_t detId, TVector3 pos, TVector3 mom, Double_t time, Double_t lenght, Double_t eloss);

  ClassDef(MpdMcord, 1) private : TClonesArray* fPointCollection;
  static MpdMcordGeo* fgGeo;
};

#endif /* MCORD_MCORD_MPDMCORD_H_ */

Zaś kod klasy przedstawiono poniżej, omówmy jego wybrane elementy.

/*
 * MpdMcord.cxx
 *
 *  Created on: 20 maj 2019
 *      Author: Daniel Wielanek
 *      E-mail: This email address is being protected from spambots. You need JavaScript enabled to view it.
 *      Warsaw University of Technology, Faculty of Physics
 */
#include "MpdMcord.h"

#include "FairGeoInterface.h"  // for FairGeoInterface
#include "FairGeoLoader.h"     // for FairGeoLoader
#include "FairGeoNode.h"       // for FairGeoNode
#include "FairGeoVolume.h"     // for FairGeoVolume
#include "FairLogger.h"        // for logging
#include "FairRootManager.h"   // for FairRootManager
#include "FairRun.h"           // for FairRun
#include "FairRuntimeDb.h"     // for FairRuntimeDb
#include "FairVolume.h"        // for FairVolume
#include "MpdDetectorList.h"   // for DetectorId::kTutDet
#include "MpdMcordGeo.h"       // for FairTutorialDet1Geo
#include "MpdMcordGeo.h"
#include "MpdMcordGeoPar.h"  // for FairTutorialDet1GeoPar
#include "MpdMcordGeoPar.h"
#include "MpdMcordPoint.h"
#include "MpdMcordPoint.h"  // for FairTutorialDet1Point
#include "MpdMcordPoint.h"
#include "MpdStack.h"
#include "TLorentzVector.h"
#include "TVirtualMC.h"

MpdMcordGeo* MpdMcord::fgGeo = NULL;


MpdMcord::MpdMcord() : MpdMcord("mcord", kTRUE) {}

MpdMcord::MpdMcord(const char* name, Bool_t active) :
  FairDetector(name, active), fTrackID(0), fVolumeID(0), fPos(0, 0, 0, 0), fMom(0, 0, 0, 0), fTime(0), fLength(0), fELoss(0) {
  fPointCollection = new TClonesArray("MpdMcordPoint");
}

MpdMcord::~MpdMcord() { delete fPointCollection; }

Jak widać głównym zadaniem konstruktora jest stworzenie tablicy z puntami MC i jej usunięcie.

Bool_t MpdMcord::ProcessHits(FairVolume* vol) {
  LOG(debug) < "In Mcord ::ProcessHits";
  // Set parameters at entrance of volume. Reset ELoss.
  if (TVirtualMC::GetMC()->IsTrackEntering()) {
    fELoss  = 0.;
    fTime   = TVirtualMC::GetMC()->TrackTime() * 1.0e09;
    fLength = TVirtualMC::GetMC()->TrackLength();
    TVirtualMC::GetMC()->TrackPosition(fPos);
    TVirtualMC::GetMC()->TrackMomentum(fMom);
  }
  // Sum energy loss for all steps in the active volume
  fELoss += TVirtualMC::GetMC()->Edep();

  // Create FairTutorialDet1Point at exit of active volume
  if (TVirtualMC::GetMC()->IsTrackExiting() || TVirtualMC::GetMC()->IsTrackStop() || TVirtualMC::GetMC()->IsTrackDisappeared()) {
    fTrackID  = TVirtualMC::GetMC()->GetStack()-GetCurrentTrackNumber();
    fVolumeID = vol->getMCid();
    if (fELoss == 0.) { return kFALSE; }
    AddHit(fTrackID, fVolumeID, fPos.Vect(), fMom.Vect(), fTime, fLength, fELoss);
  }

  return kTRUE;
}

ProcessHits to tak naprawdę "mięsko" całej klasy, to tu na podstawie informacji z Geant produkujemy punkty MC. Jest ona wywoływana za każdym razem jak cząstka wleci do detektora, rozpadnie się w nim czy będzie oddziaływać. Tutaj:

część TVirtualMC::GetMC()->IsTrackEntering() sprawdza czy cząstka pierwszy raz wlatuje do naszego detektora, jeśli tak resetujemy straty energii, ustawiamy czas wlotu cząstki, długość trajektorii oraz aktualną pozycję i pęd cząstki
fELoss += TVirtualMC::GetMC()->Edep(); będzie wywoływana gdy nasza cząstka kolejno będzie sobie przechodzić przez detektor i deponować w nim energię
jeśli jednak cząstka wyleci z detektora, zatrzyma się lub "zniknie" wtedy tworzymy dla niej punkt MC, no chyba że jej straty energetyczne są zerowe - wtedy hitu nie zapisujemy

W powyższym przykładzie kod działa tak, że zapis dokonuje się dopiero "na końcu" życia cząstki w danym detektorze, zaś zapisywany jest pęd i położenie jedynie przy wlocie do detektora.

Można jednak łatwo zmodyfikować kod tak aby zapisywał wszystkie (lub wybrane kroki symulacji) zmieniając warunek na wywołanie metody AddHit. Na koniec eventu czyścim kolekcję z punktami MC:

void MpdMcord::EndOfEvent() { fPointCollection->Delete(); }

Nie zapomnijmy jednak o zarejestrowaniu punktów MC w drzewie z danymi. W przypadku użycia wielowątkowego Geant (gMC->IsMt()==true) zapis wygląda trochę inaczej niż w zwykłym trybie.

void MpdMcord::Register() {
  if (!gMC->IsMT()) {
    FairRootManager::Instance()->Register("McordPoint", "Mcord", fPointCollection, kTRUE);
  } else {
    FairRootManager::Instance()->RegisterAny("McordPoint", fPointCollection, kTRUE);
  }
}

Istotna jest metoda konstruująca geometrię detektora:

void MpdMcord::ConstructGeometry() {
  TString fileName = GetGeometryFileName();

  if (fileName.EndsWith(".root")) {
    LOG(info) < "Constructing MCORD geometry from ROOT file" << fileName.Data();
    ConstructRootGeometry();
    return;
  }

  LOG(info) < "Constructing MCORD geometry from ASCII file  " << fileName;

  FairGeoLoader* loader          = FairGeoLoader::Instance();
  FairGeoInterface* GeoInterface = loader->getGeoInterface();
  MpdMcordGeo* MGeo              = new MpdMcordGeo();
  MGeo->setGeomFile(GetGeometryFileName());
  GeoInterface->addGeoModule(MGeo);
  Bool_t rc = GeoInterface->readSet(MGeo);
  if (rc) { MGeo->create(loader-getGeoBuilder()); }

  TList* volList = MGeo->getListOfVolumes();
  // store geo parameter
  FairRun* fRun         = FairRun::Instance();
  FairRuntimeDb* rtdb   = FairRun::Instance()->GetRuntimeDb();
  MpdMcordGeoPar* par   = static_cast<MpdMcordGeoPar*>(rtdb-getContainer("MpdMcordGeoPar"));
  TObjArray* fSensNodes = par->GetGeoSensitiveNodes();
  TObjArray* fPassNodes = par->GetGeoPassiveNodes();

  TListIter iter(volList);

  FairGeoNode* node   = nullptr;
  FairGeoVolume* aVol = nullptr;

  while ((node = (FairGeoNode*) iter.Next())) {
    aVol = dynamic_cast<FairGeoVolume*>(node);
    if (node->isSensitive())
      fSensNodes->AddLast(aVol);
    else
      fPassNodes->AddLast(aVol);
  }

  par->setChanged();
  par->setInputVersion(FairRun::Instance()-GetRunId(), 1);
  ProcessNodes(volList);
}

Kolejna ważna metoda zwraca true jeśli dany element detektora jest "aktywny". Zwykle używane jest to aby odróżnić część detektora "aktywną" tj. zbierającą dane od pozostałych elementów (np. konstrukcji podtrzymującej część aktywną, osłon).

Bool_t MpdMcord::CheckIfSensitive(std::string name) {
  TString tsname = name;
  if (tsname.BeginsWith("md01scintVol")) {
    return kTRUE;
  } else {
    return kFALSE;
  }
}

No i na koniec metoda dodająca hit. Zwrócić tu należy uwagę nie tylko na dodawanie hitu do tablicy fPointCollection ale również "informowanie" stosu z cząstkami MC (tutaj MpdStack), że nasza cząstka zostawiła jakiś sygnał w MCORD (stack->AddPoint). Jest to o tyle istotne, że zwykle cząstki które nie posiadają tej flagi są ignorowane (aby nie zajmować się zapisem cząstek MC które nie oddziaływały z detektorem).

MpdMcordPoint*
MpdMcord::AddHit(Int_t trackID, Int_t detId, TVector3 pos, TVector3 mom, Double_t time, Double_t length, Double_t eloss) {
  TClonesArray& clref = *fPointCollection;
  Int_t size          = clref.GetEntriesFast();
  MpdStack* stack     = static_cast<MpdStack*>(TVirtualMC::GetMC()->GetStack());
  stack->AddPoint(kMCORD);
  return new (clref[size]) MpdMcordPoint(trackID, detId, pos, mom, time, length, eloss);
}

Ważna informacja - zwykle we frameworkach bazujących na FairROOT do prawidłowego dodania detektora należy również zmodyfikować stos MC (w MPD to MpdStack a w CBM CbmStack) - np. w przypadku CBM natomiast AddPoint przyjmuje "enum class" który trzeba rozszerzyć o nasz detektor (ECbmModuleId). W przypadku MpdStack trzeba natomiast uważać na "filtrowanie"

void MpdMCTrack::SetNPoints(Int_t iDet, Int_t nPoints) {

  //kSTS, kTPC, kTOF, kETOF, kFD, kECT, kECAL, kNDET, kCPC, kBBC, kZDC, kFSA, kBMD,KMCORD

  Int_t mpd_nPoints = (nPoints>0)*(1 < iDet);

  if (( iDet <= kMCORD  ) && ( iDet >= kSTS  )) {
    fNPoints = ( fNPoints & ( ~ (1 < iDet) ) )  |  mpd_nPoints;
  }
  else cout < "-E- FairMCTrack::SetNPoints: Unknown detector ID "
	    < iDet << endl;

}

Z powyższego kodu widać że każdy detektor o iDet spoza zakresu kMCORD-kSTS jest odrzucany, warto dodać też nasz detektor do DetectorIdMPD (mcstack/MpdDetectorList.h).

Jeśli mamy już nase klasy to brakuje nam już tylko opisu detektora i zmuszenia MpdROOT do wstawienia naszego detektora do symulacji.

Menu

Dodawanie własnego detektora