SetMCTruth.cpp

/*
 * SetMCTruth.cpp
 *
 *  Created on: 16 Mar 2017
 *      Author: Valerio Formato
 */

/*
 * This function reads MT truth info about the primary particle and eventually
 * particle hits if they are present. Particle hits are not displayed by default
 * due to performance issues with lots of tracks in the display.
 */

#include "TEveBrowser.h"
#include "TEveScene.h"
#include "TEveTrackPropagator.h"
#include "TEveTrans.h"
#include "TGButton.h"
#include "TGTab.h"
#include "TGeoBBox.h"
#include "TGeoNode.h"
#include "TParticlePDG.h"
#include "TSystem.h"

#include "utils/GGSSmartLog.h"

#include "application/Application.h"
#include "application/event/InteractionTree.h"
#include "application/gui/MainWindow.h"

void MainWindow::SetMCTruth(GGSTMCTruthReader *reader) {
  static const std::string routineName("MainWindow::SetMCTruth");

  int nPrim = reader->GetNPrimaries();
  COUT(DEBUG) << nPrim << " Primary particles " << ENDL;
  if (nPrim == 0)
    return;

  TEveTrackList *trackList = new TEveTrackList("Primary Particles");

  TEveTrackPropagator *trkProp = trackList->GetPropagator();
  trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);

  TString primString = "";

  GGSTParticle *thisPrim = NULL;
  for (int iPrim = 0; iPrim < nPrim; iPrim++) {
    thisPrim = reader->GetPrimaryParticle(iPrim);

    COUT(DEBUG) << thisPrim << "  " << thisPrim->PDGCode << "(" << thisPrim->trackID << ")" << ENDL;

    TParticle *tempPart = new TParticle();
    tempPart->SetPdgCode(thisPrim->PDGCode);
    tempPart->SetStatusCode(1);
    tempPart->SetMother(0, 0);
    COUT(DEBUG) << tempPart->GetMass() << "  " << tempPart->GetNDaughters() << ENDL;

    TLorentzVector lpos(thisPrim->pos[0], thisPrim->pos[1], thisPrim->pos[2], thisPrim->time);
    TLorentzVector lmom;
    lmom.SetXYZM(thisPrim->mom[0], thisPrim->mom[1], thisPrim->mom[2], tempPart->GetMass());

    tempPart->SetProductionVertex(lpos);
    tempPart->SetMomentum(lmom);

    primString += "Particle: ";
    primString += tempPart->GetName();
    primString += Form("\nE_{0}: %5.3f GeV", tempPart->Energy() - tempPart->GetMass());
    primString += Form("\ntheta = %5.3f deg -  phi = %5.3f deg", 180 * (1 - tempPart->Theta() / TMath::Pi()),
                       180 * tempPart->Phi() / TMath::Pi());

    TEveTrack *thisTrack = new TEveTrack(tempPart, thisPrim->trackID, trkProp);
    thisTrack->SetLineWidth(2);
    if (tempPart->GetPDG()->Charge() > 0)
      thisTrack->SetMainColor(kCyan);
    else if (tempPart->GetPDG()->Charge() < 0)
      thisTrack->SetMainColor(kRed);
    else
      thisTrack->SetMainColor(kYellow);

    thisTrack->IncDenyDestroy();
    thisTrack->SetName(Form("Prim Track %d", thisPrim->trackID));
    thisTrack->SetStdTitle();

    trackList->AddElement(thisTrack);
  }

  trackList->MakeTracks();
  _MCTruthShortInfo->SetText(primString);

  CrossPTree();

  _MCTruthList->AddElement(trackList);
}

void MainWindow::CrossPTree() {
  static const std::string routineName("MainWindow::CrossPTree");

  EDApplication *thisApp = (EDApplication *)gApplication;

  InteractionTree<GGSTPartHit *> *pTree = thisApp->GetFileManager()->GetInteractionTree();

  TEveElementList *trackList = new TEveElementList("Particle tracks");

  TParticle *tempPart = new TParticle();

  for (std::pair<int, TreeNode<GGSTPartHit *> *> element : pTree->hashMap) {
    TEveLine *line = new TEveLine();
    if (!element.second->hits.size())
      continue;
    tempPart->SetPdgCode((int)element.second->hits[0]->particlePdg);
    if (!element.second->hits[0]->eDep)
      continue;
    for (size_t ihit = 0; ihit < element.second->hits.size(); ihit++) {
      line->SetNextPoint(element.second->hits[ihit]->entrancePoint[0], element.second->hits[ihit]->entrancePoint[1],
                         element.second->hits[ihit]->entrancePoint[2]);
      line->SetNextPoint(element.second->hits[ihit]->exitPoint[0], element.second->hits[ihit]->exitPoint[1],
                         element.second->hits[ihit]->exitPoint[2]);
    }
    line->SetName(Form("trackID %i", element.first));
    line->SetTitle(tempPart->GetName());

    if (tempPart->GetPDG()) {
      if (tempPart->GetPDG()->Charge() > 0)
        line->SetMainColor(kCyan);
      else if (tempPart->GetPDG()->Charge() < 0)
        line->SetMainColor(kRed);
      else
        line->SetMainColor(kYellow);
    } else {
      line->SetMainColor(kGray);
    }

    trackList->AddElement(line);
    // break;
  }
  // trackList->MakeTracks();
  trackList->SetRnrState(0); // Don't render
  _MCTruthList->AddElement(trackList);
  _MCPartHitList = trackList;
}