// $Id: run_g4.C,v 1.0 2007/06/01 
//
// Macro for running G4PAM_VMC with Geant4. 

void PamG4Libs_granular();

#include <PamLevel2.h>
#include <PamG4Libs.h>
#include <PamVMCApplication.h>
#include <G4ApplicationState.hh>
#include <TG4ApplicationState.h>
#include <TGeant4.h>
#include <PrimaryInfo.h>


void run_sim_g4(const TString& configMacro = "g4Config.C", const TString& primfile = "prim.root") {

  cout<<"PRIMFILE: "<<primfile<<endl;

  TFile rfin(primfile,"read");
  TTree * primtree = (TTree*)rfin.Get("primaries");
  if(!primtree) cout<<"ERROR: tree 'primaries' with primaries kinematics not found!"<<endl;
  PrimaryInfo *p = 0;
  primtree->SetBranchAddress("PRIM", &p);
  
  PamG4Libs();
    


  // Load this example libraries

  TString PAM_VMC=gSystem->Getenv("PAM_VMC");
  TString PLAT=gSystem->Getenv("PLATFORM");
  
  // MC application
  TString outname=".";
  outname+="/";
  outname+=primfile;
  outname.Replace(0,outname.Index("prim_",0),"");
  outname.Replace(outname.Index(".root",0),5,"_sim");
  
  cout<<"OUTNAME="<<outname<<endl;

  PamVMCApplication* appl 
    = new PamVMCApplication("PAMG4_VMC", "PAMELA GEANT4 VMC application", outname);

    
  // Initialize MC
  TString G4ConfMacro=PAM_VMC+"/config/"+configMacro;
  appl->InitMC(G4ConfMacro);
    
  cout<<"MACRO READED"<<endl;
    
  

  ((TGeant4*)gMC)->SetMaxNStep(50000); 

  Int_t nevents = primtree->GetEntries();
  ;
  //MAIN LOOP
  Double_t x,y,z,theta,phi,momentum;
  Int_t pdg;
  for(Int_t i = 0; i<nevents; i++){
      rfin.cd();
      primtree->GetEntry(i);
      x = p->X0;
      y = p->Y0;
      z = p->Z0;
      theta = p->Theta;
      phi = p->Phi;
      pdg = p->PDG;
      momentum = p->P0;
      
      cout<<"EVENT: "<<i<<endl;
      cout<<"PDG:"<<pdg<<endl
	  <<"X0: "<<x<<" (cm)"<<endl
	  <<"Y0: "<<y<<" (cm)"<<endl
	  <<"Z0: "<<z<<" (cm)"<<endl
	  <<"Theta:"<<theta*180./TMath::Pi()<<" (deg)"<<endl
	  <<"Phi:"<<phi*180./TMath::Pi()<<" (deg)"<<endl
	  <<"Momentum: "<<momentum<<" (GV)"<<endl;
      
      appl->GetPrimaryGenerator()->SetParticle(pdg);
      appl->GetPrimaryGenerator()->SetDirection(theta,phi);
      appl->GetPrimaryGenerator()->SetMomentum(momentum);
      appl->GetPrimaryGenerator()->SetPosition(x,y,z);

      appl->RunMC(1);

    }
    appl->FinishRun();
    //  delete appl;
}