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

void run_fluka_GP(const TString& configMacro = "flukaConfig.C") 
{
  // Load basic libraries
  gROOT->LoadMacro("./macro/basiclibs.C");
  basiclibs();

  // Load PAMELA libs
  TString PAMLIB=gSystem->Getenv("PAM_LIB");
  gSystem->Load(PAMLIB+"/libyoda.so");
  gSystem->Load(PAMLIB+"/libDarthVader.so");
  gSystem->Load(PAMLIB+"/libPamLevel2.so");


  // Load TFluka libraries
  gSystem->Load("/home/nikolas/francesco/vmc/fluka_vmc/lib/tgt_linux/libfluka");  
  // Load this example libraries
  TString G4WD=gSystem->Getenv("G4WORKDIR");
  TString PLAT=gSystem->Getenv("PLATFORM");
  gSystem->Load(G4WD+"/lib/tgt_"+PLAT+"/libPamVMC_fc");

  // MC application
  PamVMCApplication* appl 
    = new PamVMCApplication("PAMFLUKA_VMC", "PAMELA GEANT4 VMC application");
    

    appl->SetVerboseLevel(1);  
 
  
  // Initialize MC
  
  appl->InitMC(configMacro);
  // appl->SetVerboseLevel(1);  
  
  // Run MC
  

  //START READING INPUT FILE

  Int_t Ipa;
  Float_t X0, Y0, Z0, Theta, Phi, P0;
  TTree * h20 = 0;

  TString Filename = "pbar_r3_1-2.Beam.root";

  TString HOME=gSystem->Getenv("HOME");

  TSystemDirectory *readdir=new TSystemDirectory("readdir",HOME+"/RES");
    
  if ((TList*)(readdir->GetListOfFiles())->FindObject(Filename) ){
    cout<< Filename<<" found..."<<endl;
    TFile* inpf = new TFile(HOME+"/RES/"+Filename);
    h20 = (TTree*)inpf->Get("h20");
    if (h20->GetBranch("Ipa")){
      h20->SetBranchAddress("Ipa",&Ipa);
      h20->SetBranchStatus("Ipa",1);
    }
    if (h20->GetBranch("X0")){
      h20->SetBranchAddress("X0",&X0);
      h20->SetBranchStatus("X0",1);
    }
    if (h20->GetBranch("Y0")){
      h20->SetBranchAddress("Y0",&Y0);
      h20->SetBranchStatus("Y0",1);
    }
    if (h20->GetBranch("Z0")){
      h20->SetBranchAddress("Z0",&Z0);
      h20->SetBranchStatus("Z0",1);
    }
    if (h20->GetBranch("P0")){
      h20->SetBranchAddress("P0",&P0);
      h20->SetBranchStatus("P0",1);
    }
    if (h20->GetBranch("Theta")){
      h20->SetBranchAddress("Theta",&Theta);
      h20->SetBranchStatus("Theta",1);
    }
    if (h20->GetBranch("Phi")){
      h20->SetBranchAddress("Phi",&Phi);
      h20->SetBranchStatus("Phi",1);
    }
  }
  Int_t nevents = (Int_t)h20->GetEntries();

  //  nevents = 15;
  cout<<"PROCESS "<<nevents<<" EVENTS"<<endl;

  Double_t PX, PY, PZ;
  Int_t runNo = 0;
  for (Int_t i=0; i<nevents; i++){
    h20->GetEntry(i);
    cout<<"Ipa="<<Ipa<<endl;
    appl->GetPrimaryGenerator()->SetParticle(kProtonBar);
    appl->GetPrimaryGenerator()->SetMomentum(P0);
    cout<<"EKIN "<<appl->GetPrimaryGenerator()->GetKinEnergy()<<" GeV"<<endl;
    //cout<<"applP0="<<appl->GetPrimaryGenerator()->GetMomentum()<<" GV"<<endl;
    cout<<"applRIG "<<appl->GetPrimaryGenerator()->GetRigidity()<<" GV"<<endl
    PX=P0*sin(Theta)*cos(Phi);
    PY=P0*sin(Theta)*sin(Phi);
    PZ=-P0*cos(Theta);
    cout<<"(PX0,PY0,PZ0)="<<PX<<","<<PY<<","<<PZ<<endl;
    appl->GetPrimaryGenerator()->SetDirection(Theta,Phi);
    cout<<"(X0,Y0,Z0)"<<X0<<","<<Y0<<","<<Z0<<endl;
    appl->GetPrimaryGenerator()->SetPosition(X0,Y0,Z0);
    
    cout<<">StartRUN: "<<runNo<<endl;
    runNo++;

    appl->RunMC(1);

  }




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