Analiza danych cz. IV

Analiza tracków

Nadszedł czas na napisanie trochę bardziej zaawansowanego programu do analizy. A w zasadzie całej klasy będącej taskiem. Oto założenia naszego programu:

• analiza liczy rozkład pT dla danego rodzaju cząstek, każda kolekcja cząstek jest traktowana jak osobny typ cząstek
• użytkownik może analizować dowolną liczbę rodzajów cząstek - niech ta liczba będzie równa N
• możliwe, że użytkownik doda kilka (M) kolekcji cząstek, wtedy dla każdego zderzenia liczymy M-kolekcji cząstek, musimy więc zapisać NxM histogramów z rozkładem pT

Nasza analiza bazuje na właściwościach cząstek, więc naszą bazową klasą powinna być Hal::TrackAna. Napiszemy kod na dwa sposoby (ale wynik analizy będą takie same).

• opcja pierwsza - każda kolekcja cząstek jest replikowana N-razy, mamy więc efektywnie NxM kolekcji cząstek i każda ma swój histogram
• opcja druga - nie replikujemy kolekcji

Kiedy użyć pierwszej opcji a kiedy drugiej? Poniżej połączenia kolekcji w opcji pierwszej:

oraz w opcji drugiej:

Na czym polega różnica? na liczbie kolekcji i połączeniach między nimi. Załóżmy że kolekcja zderzeń nr 0 to zderzenia centralne, a kolekcja 1 to zderzenia peryferyjne, zaś w przypadku kolekcji mamy kolejno pion, kaony i protony. Jak wygląda przepływ danych?

• w opcji pierwsze eventy z kolekcji 0 są źródłem cząstek dla kolekcji 0 1 i 2 a eventy z kolekcji 1 dla cząstek z kolekcji 3,4,5
• w opcji drugiej każda kolekcja eventów generuje cząstki dla każdej kolekcji

W praktyce oznacza to że jeśli np. kolekcja 0 to zderzenia centralne, a kolekcja 1 to zderzenia peryferyjne, zaś kolekcje cząstek to odpowiednio piony,kaony i protony to nasze histogramy będą zawierały następujące dane:

Zacznę od opisu klasy, gdzie replikujemy cięcia cząstkowe (opcja pierwsza). Na samym końcu strony jest pełen kod więc skupię się tylko na wybranych metodach. Zacznę od opcji pierwszej z replikacja kolekcji cząstek.

Funkcja Init:

• wywołujemy metodę Init klasy nadrzędnej, generalnie inicjalizacja tasku z analizą to dużo rzeczy więc polecam robić inicjalizację klasy nadrzędnej a dopiero potem dodanych przez siebie rzeczy
• sprawdzamy czy mamy co najmniej jedną kolekcję cząstek i zderzeń, jeśli nie, to zwracamy kERROR i task nie będzie uruchomiony
• tworzymy HistogramAxisConf - to klasa opisująca osie histogramu
• tworzymy HistogramManager - to klasa będąca tablicą obiektów-histogramów

  virtual Hal::Task::EInitFlag Init() {
    auto status    = Hal::TrackAna::Init();
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    if (eventCol == 0) return Hal::Task::EInitFlag::kERROR;
    if (trackCol == 0) return Hal::Task::EInitFlag::kERROR;
    std::vector<Hal::Histogramaxisconf> axisconf;
    Hal::HistogramAxisConf xAxis("pT [GeV/c]", fBins, fMin, fMax);
    Hal::HistogramAxisConf yAxis("N [tracks]", 0, 0, 0);
    axisconf.push_back(xAxis);
    axisconf.push_back(yAxis);
    fHistograms = new Hal::HistogramManager_1_1D<TH1D>();
    fHistograms->Init(trackCol, axisconf, kTRUE);
    int mod   = trackCol / eventCol;
    int count = 0;
    for (int i = 0; i < eventCol; i++) {
      for (int j = 0; j < mod; j++) {
        fHistograms->At(count++)->SetTitle(Form("%i %i", i, j));
      }
    }
    return status;
  }  

Kolejna funkcja to ProcessTrack, jest ona wywołna dla każdej cząstki, która przeszła cięcia. Tutaj wyliczamy pęd poprzeczny i wypełniamy nim histogram. Indeks histogramu to także numer aktualnie przeliczanej kolekcji cząstek. 

protected:
  virtual void ProcessTrack() {
    double pt = fCurrentTrack->GetMomentum().Pt();
    fHistograms->Fill(fCurrentTrackCollectionID, pt);
  }

Kolejna metoda to LinkConnections, ona definuje jak łączą się nasze kolekcje, ponieważ kolekcji cząstek mamy NxM, to każdą kolekcję zderzeń łączymy z M kolekcją cząstek. Zakomentowana linia robi to samo co pętla for iEvent - replikuje kolekcje cząstkowe N-razy.

  virtual void LinkCollections() {//też protected
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();

    // fCutContainer->LinkCollections(Hal::ECutUpdate::kEvent, Hal::ECutUpdate::kTrack,
    // Hal::CutContainer::ELinkPolicy::kReplicateLast);
    Int_t mod = trackCol / eventCol;
    for (int iEvent = 0; iEvent < eventCol; iEvent++) {
      for (int iTrack = 0; iTrack < mod; iTrack++) {
        fCutContainer->LinkCollections(Hal::ECutUpdate::kEvent, iEvent, Hal::ECutUpdate::kTrack, iTrack + mod * iEvent);
      }
    }
  }  

CheckCutContainersCollections sprawdza czy liczba kolekcji jest ok. Tutaj replikujemy nasze kolekcje cząstek N razy (liczba kolekcji zderzeń). Jest ona wywoływana przed LinkCollections.

  virtual void CheckCutContainerCollections() {//protected
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    for (int iEvent = 1; iEvent < eventCol; iEvent++) {
      for (int iTrack = 0; iTrack < trackCol; iTrack++) {
        fCutContainer->ReplicateCollection(Hal::ECutUpdate::kTrack, iTrack, trackCol * iEvent + iTrack);
      }
    }
  }

Ostatnia metoda to Report. Tutaj tylko dodajemy obiekty za analizy do paczki danych. Paczka musi mieć określoną strukturę (inaczej wszystkie raporty czy dostęp do histogramów przez Hal:AnaFile mogą nie działać prawidłowo), dlatego nie tworzymy jej od zera, ale używamy paczki stworzonej przez klasę nadrzędną i tylko dorzucamy nasze histogramy. Paczka jest tym co fizycznie będzie zapisane w pliku wynikowym.

  virtual Hal::Package* Report() const {//protected
    auto report = Hal::TrackAna::Report();
    for (int i = 0; i < fHistograms->GetSize(); i++) {
      report->AddObject(fHistograms->At(i)->Clone());
    }
    return report;
  }  

W przypadku opcji drugiej niemal wszystko wygląda tak samo z kilkoma drobnymi różnicami, wynikają one z tego że inaczej trzeba przeliczać numery kolekcji i tak np. zmieniony Init wygląda tak:

    virtual Hal::Task::EInitFlag Init() {
    auto status    = Hal::TrackAna::Init();
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    if (eventCol == 0) return Hal::Task::EInitFlag::kERROR;
    if (trackCol == 0) return Hal::Task::EInitFlag::kERROR;
    std::vector<Hal::Histogramaxisconf> axisconf;
    Hal::HistogramAxisConf xAxis("pT [GeV/c]", fBins, fMin, fMax);
    Hal::HistogramAxisConf yAxis("N [tracks]", 0, 0, 0);
    axisconf.push_back(xAxis);
    axisconf.push_back(yAxis);
    fHistograms = new Hal::HistogramManager_1_1D<TH1D>();
    fHistograms->Init(trackCol*eventCol, axisconf, kTRUE);
    int mod   = trackCol / eventCol;
    int count = 0;
    for (int i = 0; i < eventCol; i++) {
      for (int j = 0; j < mod; j++) {
        fHistograms->At(count++)->SetTitle(Form("%i %i", i, j));
      }
    }
    return status;
  }

W metodzie ProcessTrack również nieco inaczej trzeba przeliczać numer wypełnianego histogramu, bo mamy więcej histogramów (MxN) niż kolekcji cząstek (M). Jeśli wypełnimy histogram tak jak w poprzednim kodzie - używają jako indeksu fCurrentTrackCollectioID nie będziemy w stanie odróżnić cząstek z poszczególnych kolekcji zderzeń bo dla takiej samej kolekcji cząstek możemy mieć dane z różnych kolekcji zderzeń.

    virtual void ProcessTrack() {
    double pt = fCurrentTrack->GetMomentum().Pt();
    fHistograms->Fill(fTrackCollectionsNo * fCurrentEventCollectionID + fCurrentTrackCollectionID, pt);
  }

Linkowanie również trzeba zrobić inaczej, w tym przypadku użyłem innej metody, która automatycznie linkuje każdą kolekcję zderzeń z każdą kolekcję cząstek. 

    virtual void LinkCollections() {
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    Int_t mod      = trackCol / eventCol;
    fCutContainer->LinkCollections(Hal::ECutUpdate::kEvent, Hal::ECutUpdate::kTrack, Hal::CutContainer::ELinkPolicy::kAnyToAny);
  }
  

Poniżej kompletny kod, dwa różne podejścia są przedstawione jako dwie różne klasy MyAna1 i MyAna2:

 class MyAna : public Hal::TrackAna {
  Hal::HistogramManager_1_1D<TH1D>* fHistograms = {nullptr};
  Double_t fMin = {0.0}, fMax = {1.0};
  Int_t fBins = {100};

public:
  MyAna() {};
  void ConfigAxis(Int_t bins, Double_t lo, Double_t hi) {
    fBins = bins;
    fMin  = lo;
    fMax  = hi;
  }

  virtual ~MyAna() {
    if (fHistograms) delete fHistograms;
  }
  virtual Hal::Task::EInitFlag Init() {
    auto status    = Hal::TrackAna::Init();
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    if (eventCol == 0) return Hal::Task::EInitFlag::kERROR;
    if (trackCol == 0) return Hal::Task::EInitFlag::kERROR;
    std::vector<Hal::HistogramAxisConf> axisconf;
    Hal::HistogramAxisConf xAxis("pT [GeV/c]", fBins, fMin, fMax);
    Hal::HistogramAxisConf yAxis("N [tracks]", 0, 0, 0);
    axisconf.push_back(xAxis);
    axisconf.push_back(yAxis);
    fHistograms = new Hal::HistogramManager_1_1D<TH1D>();
    fHistograms->Init(trackCol, axisconf, kTRUE);
    int mod   = trackCol / eventCol;
    int count = 0;
    for (int i = 0; i < eventCol; i++) {
      for (int j = 0; j < mod; j++) {
        fHistograms->At(count++)->SetTitle(Form("%i %i", i, j));
      }
    }
    return status;
  }

protected:
  virtual void ProcessTrack() {
    double pt = fCurrentTrack->GetMomentum().Pt();
    fHistograms->Fill(fCurrentTrackCollectionID, pt);
  }

  virtual void LinkCollections() {
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();

    // fCutContainer->LinkCollections(Hal::ECutUpdate::kEvent, Hal::ECutUpdate::kTrack,
    // Hal::CutContainer::ELinkPolicy::kReplicateLast);
    Int_t mod = trackCol / eventCol;
    for (int iEvent = 0; iEvent < eventCol; iEvent++) {
      for (int iTrack = 0; iTrack < mod; iTrack++) {
        fCutContainer->LinkCollections(Hal::ECutUpdate::kEvent, iEvent, Hal::ECutUpdate::kTrack, iTrack + mod * iEvent);
      }
    }
  }
  virtual Hal::Package* Report() const {
    auto report = Hal::TrackAna::Report();
    for (int i = 0; i < fHistograms->GetSize(); i++) {
      report->AddObject(fHistograms->At(i)->Clone());
    }
    return report;
  }
  virtual void CheckCutContainerCollections() {
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    for (int iEvent = 1; iEvent < eventCol; iEvent++) {
      for (int iTrack = 0; iTrack < trackCol; iTrack++) {
        fCutContainer->ReplicateCollection(Hal::ECutUpdate::kTrack, iTrack, trackCol * iEvent + iTrack);
      }
    }
  }
};

class MyAna2 : public Hal::TrackAna {
  Hal::HistogramManager_1_1D<TH1D>* fHistograms = {nullptr};
  Double_t fMin = {0.0}, fMax = {1.0};
  Int_t fBins = {100};

public:
  MyAna2() {};
  void ConfigAxis(Int_t bins, Double_t lo, Double_t hi) {
    fBins = bins;
    fMin  = lo;
    fMax  = hi;
  }

  virtual ~MyAna2() {
    if (fHistograms) delete fHistograms;
  }
  virtual Hal::Task::EInitFlag Init() {
    auto status    = Hal::TrackAna::Init();
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    if (eventCol == 0) return Hal::Task::EInitFlag::kERROR;
    if (trackCol == 0) return Hal::Task::EInitFlag::kERROR;
    std::vector<Hal::HistogramAxisConf> axisconf;
    Hal::HistogramAxisConf xAxis("pT [GeV/c]", fBins, fMin, fMax);
    Hal::HistogramAxisConf yAxis("N [tracks]", 0, 0, 0);
    axisconf.push_back(xAxis);
    axisconf.push_back(yAxis);
    fHistograms = new Hal::HistogramManager_1_1D<TH1D>();
    fHistograms->Init(eventCol * trackCol, axisconf, kTRUE);
    int mod   = trackCol / eventCol;
    int count = 0;
    for (int i = 0; i < eventCol; i++) {
      for (int j = 0; j < mod; j++) {
        fHistograms->At(count++)->SetTitle(Form("%i %i", i, j));
      }
    }
    return status;
  }

protected:
  virtual void ProcessTrack() {
    double pt = fCurrentTrack->GetMomentum().Pt();
    fHistograms->Fill(fTrackCollectionsNo * fCurrentEventCollectionID + fCurrentTrackCollectionID, pt);
  }

  virtual void LinkCollections() {
    Int_t eventCol = fCutContainer->GetEventCollectionsNo();
    Int_t trackCol = fCutContainer->GetTrackCollectionsNo();
    Int_t mod      = trackCol / eventCol;
    fCutContainer->LinkCollections(Hal::ECutUpdate::kEvent, Hal::ECutUpdate::kTrack, Hal::CutContainer::ELinkPolicy::kAnyToAny);
  }
  virtual Hal::Package* Report() const {
    auto report = Hal::TrackAna::Report();
    for (int i = 0; i < fHistograms->GetSize(); i++) {
      report->AddObject(fHistograms->At(i)->Clone());
    }
    return report;
  }
};


void AddCuts(Hal::TrackAna* ana) {
  Hal::EventPhiCut phi;
  phi.SetMinMax(0, TMath::Pi());
  ana->AddCut(phi, "{0}");
  phi.SetMinMax(-TMath::Pi(), TMath::Pi());
  ana->AddCut(phi, "{1}");
  Hal::TrackPtCut pT;
  pT.SetMinMax(0, 0.25);
  ana->AddCut(pT, "{0}");
  pT.SetMinMax(0.25, 0.5);
  ana->AddCut(pT, "{1}");
  pT.SetMinMax(0.5, 1.0);
  ana->AddCut(pT, "{2}");
}

void track_ana() {
  auto run       = new Hal::AnalysisManager();
  auto source    = new HalOTF::Source(10000);
  auto generator = new HalOTF::EventGenerator();
  TH2D h("spec", "spec", 1, -1, 1, 1, 0, 1);
  h.SetBinContent(1, 1, 1);
  generator->SetSpiecies(h, Hal::Const::PionPlusPID());
  generator->SetFixMult(100);

  source->AddEventGenerator(generator);
  source->Register();
  run->SetOutput("/opt/temp/glauber.root");
  run->SetSource(source);
  auto reader = new HalOTF::Reader();
  run->AddTask(reader);

  MyAna* ana1 = new MyAna();
  AddCuts(ana1);
  MyAna2* ana2 = new MyAna2();
  AddCuts(ana2);
  run->AddTask(ana1);
  run->AddTask(ana2);

  run->Init();
  run->Run();
}