// $Id: run_g4.C,v 1.0 2007/06/01 
//
// Macro for running G4PAM_VMC with Geant4. 
void run_g4_GP(const TString& configMacro = "g4Config.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 Geant4 libraries
    
    gROOT->LoadMacro("./macro/PamG4Libs.C");
    PamG4Libs();
    
    // 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("PAMG4_VMC", "PAMELA GEANT4 VMC application");  
     
    
    // Initialize MC
    appl->InitMC(configMacro);
    ((TGeant4*)gMC)->ProcessGeantMacro("g4config2.in");

    //START READING INPUT FILE

    Int_t Ipa;
    Float_t X0, Y0, Z0, Theta, Phi, P0;

    TTree * h20 = 0;

    TString Filename = "protons_N0_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 = 1000; 
    cout<<"PROCESS "<<nevents<<" EVENTS"<<endl;

    TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(kProton); //correct
    Double_t mass = particlePDG->Mass();
    Double_t PX, PY, PZ;
    for (Int_t i=0; i<nevents; i++){
      h20->GetEntry(i);
      cout<<"Ipa="<<Ipa<<endl;
      appl->GetPrimaryGenerator()->SetParticle(kProton);
      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);

      
            appl->RunMC(1);
     
    }

  
   // Visualization settings
    // ((TGeant4*)gMC)->ProcessGeantMacro("g4vis.in");
    // inpf.cd();
    //   delete inpf;
    // Run MC
    //  appl->RunMC(1000);
      appl->FinishRun();
      ((TGeant4*)gMC)->StartGeantUI();
      //      delete appl;

  }