ExtractStoppingPower.cc

////////////////////////////////////////////////////////
/// Extracts the stopping power of the primary track.
/// P G Jones <[email protected]>
///                                 
/// 02/01/2013 : New File                             
////////////////////////////////////////////////////////

#include <RAT/DS/Root.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/Run.hh>

#include <TGraph.h>
#include <TTree.h>
#include <TFile.h>

#include <time.h>
using namespace std;

int graphPoint = 0;

void
FillHistogram( RAT::DS::MC* rMC,
               TGraph* nCer );

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun );

void
ExtractStoppingPower( const char* inFile,
                      const char* outFile )
{
  // Load the root file first
  RAT::DS::Root* rDS;
  RAT::DS::Run* rRun;
  TTree *tree;
  LoadRootFile( inFile, &tree, &rDS, &rRun );

  time_t codeStart = time( NULL );
  TGraph* dedx = new TGraph();

  for( int iEvent = 0; iEvent < tree->GetEntries(); iEvent++ )
    {
      if( iEvent % 100 == 0 )
        cout << iEvent << " finished at " << time( NULL ) - codeStart << endl;
      tree->GetEntry( iEvent );
      FillHistogram( rDS->GetMC(), dedx );
    }

  TFile outputFile( outFile, "RECREATE" );
  outputFile.cd();
  dedx->SetName( "dedx" );
  dedx->Write();
  outputFile.Close();
}

void
FillHistogram( RAT::DS::MC* rMC,
               TGraph* dedx )
{
  for( int iTrack = 0; iTrack < rMC->GetMCTrackCount(); iTrack++ )
    {
      RAT::DS::MCTrack* rMCTrack = rMC->GetMCTrack( iTrack );
      if( rMCTrack->GetParentID() != 0 ) // Primary track has parent ID 0
        continue;
      double previousEnergy = rMCTrack->GetMCTrackStep( 0 )->GetKE();
      for( int iStep = 1; iStep < rMCTrack->GetMCTrackStepCount(); iStep++ )
  {
          RAT::DS::MCTrackStep* rMCTrackStep = rMCTrack->GetMCTrackStep( iStep );
          double dE = previousEnergy - rMCTrackStep->GetKE();
          if( dE > 0.0 && rMCTrackStep->GetLength() > 0.0 )
            {
              cout << rMCTrackStep->GetKE() << " " << dE << " " << rMCTrackStep->GetLength() << endl;
              dedx->SetPoint( graphPoint++, previousEnergy, 10.0 * dE / rMCTrackStep->GetLength() ); // Convert to cm
            }
          previousEnergy = rMCTrackStep->GetKE();
	}
    }
}

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun )
{
  TFile *file = new TFile( lpFile );
  (*tree) = (TTree*)file->Get( "T" );
  TTree *runTree = (TTree*)file->Get( "runT" );

  *rDS = new RAT::DS::Root();

  (*tree)->SetBranchAddress( "ds", &(*rDS) );

  *rRun = new RAT::DS::Run();

  runTree->SetBranchAddress( "run", &(*rRun) );
  runTree->GetEntry();
}