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 "TSystem.h"
#include "TGTab.h"
#include "TGButton.h"
#include "TGeoBBox.h"
#include "TGeoNode.h"
#include "TEveBrowser.h"
#include "TEveScene.h"
#include "TEveTrans.h"
#include "TEveTrackPropagator.h"
#include "TParticlePDG.h"

#include "utils/GGSSmartLog.h"

#include "application/gui/MainWindow.h"
#include "application/event/InteractionTree.h"
#include "application/Application.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;
}