// $Id: run_g4.C,v 1.0 2007/06/01 
//
// Macro for running G4PAM_VMC with Geant4. 
void run_g4_GP_Chain(const TString& configMacro = "../macro/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");
   
    appl->GetPrimaryGenerator()->SetNofPrimaries(1);  
     
    
    // Initialize MC
    appl->InitMC(configMacro);
    ((TGeant4*)gMC)->ProcessGeantMacro("../macro/g4config2.in");

    //START READING INPUT FILE

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

    TChain * h20 = new TChain("h20");

    //TString Filename = "ele_5-6_n2.Beam.root";

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

    TSystemDirectory *readdir=new TSystemDirectory("readdir",HOME+"/RES");

    // TList * list = readdir->GetListOfFiles();

    h20->Add(HOME+"/RES/ele_r1_5-6.Beam.root");
    h20->Add(HOME+"/RES/ele_r2_5-6.Beam.root");
    h20->Add(HOME+"/RES/ele_r3_5-6.Beam.root");
    //    h20->Add(HOME+"/RES/ele_n4_5-6.Beam.root");
    //if (list->FindObject(Filename) ){
    // cout<< Filename<<" found..."<<endl;
      //   TFile* inpf = new TFile(HOME+"/RES/"+Filename);

      //h20 = (Te*)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 = 200; 
    cout<<"PROCESS "<<nevents<<" EVENTS"<<endl;

    TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(kElectron); //correct
    Double_t mass = particlePDG->Mass();
    Double_t PX, PY, PZ;
    for (Int_t i=0; i<nevents; i++){
      h20->GetEntry(i);
      appl->GetPrimaryGenerator()->SetKinEnergy(sqrt(P0*P0+mass*mass)-mass);
      cout<<"EKIN "<<appl->GetPrimaryGenerator()->GetKinEnergy()<<" GeV"<<endl;
      //cout<<"RIG "<<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(PX,PY,PZ);
      cout<<"(X0,Y0,Z0)"<<X0<<","<<Y0<<","<<Z0<<endl;
      appl->GetPrimaryGenerator()->SetPosition(X0,Y0,Z0);

      
            appl->RunMC(1);
     
    }

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

  }