//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: G4PiMinusAbsorptionAtRestTest.cc,v 1.5 2008/12/18 13:02:40 gunter Exp $
// GEANT4 tag $Name: geant4-09-03-cand-01 $
//
// -------------------------------------------------------------------
//      GEANT 4 class file --- Copyright CERN 1998
//      CERN Geneva Switzerland
//
//
//      File name:     G4PiMinusAbsorptionAtRestTest.cc 
//
//      Author:        Maria Grazia Pia (pia@genova.infn.it), 
//                     (from Ch. Voelcker, from M. Maire)
// 
//      Creation date: 8 May 1998
//
//      Modifications: 
//      MG Pia         6 Jul 1998 Modified handling of ProcessManager
//                                to be consistent with changes in 
//                                ParticleDefinition
//
// -------------------------------------------------------------------

#include "G4ios.hh"
#include <fstream>
#include <iomanip>

#include "G4Material.hh"
#include "G4Element.hh"

#include "G4ProcessManager.hh"
#include "G4PiMinusAbsorptionAtRest.hh"
#include "G4PionMinusAbsorptionAtRest.hh"

#include "G4ParticleDefinition.hh"
#include "G4DynamicParticle.hh"
#include "G4PionMinus.hh"
#include "G4ParticleMomentum.hh"

#include "G4Box.hh"
#include "G4PVPlacement.hh"

#include "G4GRSVolume.hh"
#include "G4LogicalVolume.hh"

#include "G4ProcessManager.hh"
#include "G4ForceCondition.hh"

#include "G4Step.hh"
#include "G4Track.hh"

#include "CLHEP/Hist/TupleManager.h"
#include "CLHEP/Hist/HBookFile.h"
#include "CLHEP/Hist/Histogram.h"
#include "CLHEP/Hist/Tuple.h"

#include "G4IonTable.hh"
#include "G4Electron.hh"
#include "G4Proton.hh"
#include "G4Neutron.hh"
#include "G4ParticleTypes.hh"
#include "G4IonTable.hh"
#include "G4ParticleTable.hh"
#include "G4GenericIon.hh"

int main()
{
  G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
  G4ParticleDefinition* ion  = G4GenericIon::GenericIonDefinition();

  G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
  proton->SetCuts(0.01);
  G4ParticleDefinition* neutron  = G4Neutron::NeutronDefinition();
  neutron->SetCuts(0.01);


  G4ParticleTable* theTableOfParticles;
  theTableOfParticles = G4ParticleTable::GetParticleTable();
  G4IonTable* theTable = new G4IonTable();

  // -------------------------------------------------------------------
  // MGP ---- HBOOK initialization
  HepTupleManager* hbookManager;
  hbookManager = new HBookFile("mg.hbook", 58);
  assert (hbookManager != 0);

  // MGP ---- Book a histogram
  HepHistogram* hEKin;
  hEKin = hbookManager->histogram("Kinetic Energy", 100,0.,150.);
  assert (hEKin != 0);  

  HepHistogram* hP;
  hP = hbookManager->histogram("Momentum", 100,0.,1000.);
  assert (hP != 0);  

  HepHistogram* hNSec;
  hNSec = hbookManager->histogram("Number of secondaries", 25,0.,25.);
  assert (hNSec != 0);  

  HepHistogram* hDebug;
  hDebug = hbookManager->histogram("Debug", 100,0.,140.);
  assert (hDebug != 0);  

  // MGP ---- Book a ntuple
  HepTuple* ntuple;
  ntuple = hbookManager->ntuple("Pion absorption at rest Ntuple");
  assert (ntuple != 0);

  G4String name, symbol;
  G4double a, iz, z, density;
  G4int nEl;

  G4int nIter = 10;
  G4int imat = 3;
  G4int verboseLevel = 1;
  G4int processId = 1;
  G4int deexcitationIndex =0;

  cout << " 0) Copper      1) Lead         2) Iron         3) Carbon" << G4endl;
  cout << " 4) LArgon      5) Polystyrene  6) Tungsten     7) Oxygen" << G4endl;
  cout << " 8) Beryllium   9) Aluminium   10) Uranium     11) BGO" << G4endl;
  cout << "12) NaI        13) CsI         14) Kapton" << G4endl;

  G4cout 
    << "Enter number of absorptions [10], material [3], Verbose level [1]"
    << G4endl;
  G4cin >> nIter >> imat >> verboseLevel;

  G4cout 
    << "Enter process: 1 G4PiMinusAbsorptionAtRest, 2 G4GheishaAbsAtRest" 
    << G4endl;
  G4cin >> processId;

  if (processId < 1 || processId >2)
    {
      G4cout << "Wrong process ID, set to default" << G4endl;
      processId = 1;
    }

  G4cout << "Enter deexcitation algorithm: 0 Theo, 1 Dummy" << G4endl;
  G4cin >> deexcitationIndex;
  
  G4int hnt;
  G4cout << "Enter histo (0) or ntuple (1) " << G4endl;
  G4cin >> hnt;

  // Materials definition 

  // Materials definition 

    G4Material* Cu = new G4Material(name="Copper", density=8.96*g/cm3, nEl=1);
    G4Element* elCu = new G4Element(name="Copper", symbol="Cu", iz=29., a=63.55*g/mole);
    Cu->AddElement( elCu, 1 );
    
    G4Material* Pb = new G4Material(name="Lead", density=11.35*g/cm3, nEl=1);
    G4Element* elPb = new G4Element(name="Lead", symbol="Pb", iz=82., a=207.19*g/mole);
    Pb->AddElement( elPb, 1 );
    
    G4Material*  Fe = new G4Material(name="Iron", density=7.87*g/cm3, nEl=1);
    G4Element* elFe = new G4Element(name="Iron", symbol="Fe", iz=26., a=55.85*g/mole);
    Fe->AddElement( elFe, 1 );
    
    G4Material*  Graphite = new G4Material(name="Graphite", density=2.265*g/cm3, nEl=1);
    G4Element* elC  = new G4Element(name="Carbon", symbol="C", iz=6., a=12.0107*g/mole);
    Graphite->AddElement( elC, 1 );

    G4Material*  LAr= new G4Material(name="LArgon", density=1.393*g/cm3, nEl=1);
    G4Element* elAr  = new G4Element(name="Argon", symbol="Ar", iz=18., a=39.95*g/mole);
    LAr->AddElement( elAr, 1 );
    
    G4Material*  PS = new G4Material(name="PolyStyrene", density=1.032*g/cm3, nEl=2);
    //    G4Element* elC = new G4Element(name="Carbon", symbol="C", iz=6., a=12.01*g/mole);
    G4Element* elH = new G4Element(name="Hydrogen", symbol="H", iz=1., a=1.01*g/mole);
    PS->AddElement( elC, 8 );
    PS->AddElement( elH, 8 );
    
    G4Material*  W  = new G4Material(name="Tungsten", density=19.30*g/cm3, nEl=1);
    G4Element* elW  = new G4Element(name="Tungsten", symbol="W", iz=74., a=183.85*g/mole);
    W->AddElement( elW, 1 );
    
    // approximate numbers for O
    G4Material*  O = new G4Material(name="Oxygen", density=1.1*g/cm3, nEl=1);
    G4Element* elO  = new G4Element(name="Oxygen", symbol="O", iz=8., a=15.9994*g/mole);
    O->AddElement( elO,  1 );
    
    G4Material*  Be = new G4Material(name="Beryllium", density=1.848*g/cm3, nEl=1);
    G4Element* elBe  = new G4Element(name="Beryllium", symbol="Be", iz=4., a=9.01*g/mole);
    Be->AddElement( elBe, 1 );
    
    G4Material*  Al = new G4Material(name="Aluminium", density=2.70*g/cm3, nEl=1);
    G4Element* elAl  = new G4Element(name="Aluminium", symbol="Al", iz=13., a=26.98*g/mole);
    Al->AddElement( elAl, 1 );
    
    
    G4Material*  U = new G4Material(name="Uranium", density=18.95*g/cm3, nEl=1);
    G4Element* elU  = new G4Element(name="Uranium", symbol="U", iz=92., a=238.03*g/mole);
    U->AddElement( elU, 1 );
    
    G4Material*  BGO = new G4Material(name="BGO", density=2.15*g/cm3, nEl=3);
    G4Element* elBi = new G4Element(name="Bismuth", symbol="Bi", iz=83., a=208.98*g/mole);
    G4Element* elGe = new G4Element(name="Germanium", symbol="Ge", iz=32., a=72.59*g/mole);
    BGO->AddElement( elBi, 4 );
    BGO->AddElement( elGe, 3 );
    BGO->AddElement( elO, 12 );
    
    G4Material*  NaI = new G4Material(name="NaI", density=3.67*g/cm3, nEl=2);
    G4Element* elNa = new G4Element(name="Sodium", symbol="Na", iz=11., a=22.990*g/mole);
    G4Element* elI = new G4Element(name="Iodine", symbol="I", iz=53., a=126.904*g/mole);
    NaI->AddElement( elNa, 1 );
    NaI->AddElement( elI, 1 );
    
    G4Material*  CsI = new G4Material(name="CsI", density=4.53*g/cm3, nEl=2);
    G4Element* elCs = new G4Element(name="Cesium", symbol="Cs", iz=55., a=132.905*g/mole);
    CsI->AddElement( elCs, 1 );
    CsI->AddElement( elI, 1 );

    G4Material*  Kapton = new G4Material(name="Kapton", density=1.53*g/cm3, nEl=4); 
    // formula: private communications, see mail.
    Kapton->AddElement( elC, 22 );
    Kapton->AddElement( elH, 10 );
    Kapton->AddElement( elO,  5 );
    G4Element* elN = new G4Element(name="Nitrogen", symbol="N", iz=7., a=14.007*g/mole);
    Kapton->AddElement( elN, 2 );


    //  G4Material* Al = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3 );
    //  G4Material* Fe = new G4Material("Iron",      26., 55.85*g/mole, 7.87*g/cm3 );
    //  G4Material* Pb = new G4Material("Lead",      82., 207.19*g/mole, 11.35*g/cm3 );
    //  G4Material* Graphite = new G4Material("Graphite", 6., 12.00*g/mole, 2.265*g/cm3 );
  
    //  G4Element*   H = new G4Element ("Hydrogen", "H", 1. ,  1.01*g/mole);
    //  G4Element*   O = new G4Element ("Oxygen"  , "O", 8. , 16.00*g/mole);
    //  G4Element*   C = new G4Element ("Carbon"  , "C", 6. , 12.00*g/mole);
    //  G4Element*  Cs = new G4Element ("Cesium"  , "Cs", 55. , 132.905*g/mole);
    //  G4Element*   I = new G4Element ("Iodide"  , "I", 53. , 126.9044*g/mole);
  
    //  G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
    //  water->AddElement(H,2);
    //  water->AddElement(O,1);
    //  
    //  G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
    //  ethane->AddElement(H,6);
    //  ethane->AddElement(C,2);
  
    //  G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
    //  csi->AddElement(Cs,1);
    //  csi->AddElement(I,1);
  
  //  G4Element::DumpInfo(); 
  //  G4Material::DumpInfo();
  
  static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
  
  // Fix a material
  // G4int imat = 0;   // graphite
  
  // Geometry definitions

  G4Box* theFrame = new G4Box ("Frame",1*m, 1*m, 1*m);
  G4LogicalVolume* LogicalFrame = 
    new G4LogicalVolume (theFrame,(*theMaterialTable)(imat),"LFrame", 0,0,0);
  
  G4PVPlacement* PhysicalFrame = 
    new G4PVPlacement(0,G4ThreeVector(),
		      "PFrame",LogicalFrame,0,false,0);
  
  // The center-of-mass of the cube should be located at the origin!

  // Particle 
  G4ParticleDefinition* pionMinus = G4PionMinus::PionMinusDefinition();

  // Processes
  G4ProcessManager* thePionProcessManager = new G4ProcessManager(pionMinus);
  pionMinus->SetProcessManager(thePionProcessManager);
  
  //  thePionProcessManager->SetParticleType(pionMinus);

  G4PiMinusAbsorptionAtRest* thePiAbsorptionProcess = 
    new G4PiMinusAbsorptionAtRest("G4PiMinusAbsorptionAtRest");
  G4PionMinusAbsorptionAtRest* gheishaAbsorptionProcess =
    new G4PionMinusAbsorptionAtRest();

  switch (processId)
    {
    case 1:
      {
	thePionProcessManager->AddProcess(thePiAbsorptionProcess,0,-1,0);
	thePiAbsorptionProcess->SetDeexcitationAlgorithm(deexcitationIndex);  
	break;
      }
    case 2:
      {
	thePionProcessManager->AddProcess(gheishaAbsorptionProcess,0,-1,0);
        break;
      }
    default:
      {
	thePionProcessManager->AddProcess(thePiAbsorptionProcess,0,-1,0);
	thePiAbsorptionProcess->SetDeexcitationAlgorithm(deexcitationIndex);  
	break;
      }
    }

  G4ForceCondition* condition;
  
  // Set cut and Build CrossSection Tables 
  pionMinus->SetCuts(1.*mm);
  
  // Create one Dynamic Particle  
  G4double pionEnergy = 0.*MeV;
  G4ParticleMomentum pionDirection(0.,0.,1.);
  G4DynamicParticle aPionMinus(G4PionMinus::PionMinus(),pionDirection,pionEnergy);
  
  // Track definition (for this test ONLY!)
  G4ThreeVector aPosition(0.,0.,0.);
  G4double aTime = 0. ;
  G4Track* ptrack = new G4Track(&aPionMinus,aTime,aPosition) ;
  G4Track& aTrack = (*ptrack) ;
  
  //ptrack->SetVolume(PhysicalFrame);
  
  // Do I really need this?
  G4GRSVolume* touche = new G4GRSVolume(PhysicalFrame, NULL, aPosition);   
  ptrack->SetTouchable(touche);
  
  // Create 1 Step (for this test only)
  G4Step* Step = new G4Step();  
  G4Step& aStep = (*Step);
  Step->SetTrack(ptrack);
  
  // Check applicability
  G4ParticleDefinition* PionMinusDefinition = aPionMinus.GetDefinition();

  if (! thePiAbsorptionProcess->IsApplicable(*PionMinusDefinition)) 
    {
      G4cout 
	<< PionMinusDefinition->GetParticleName() << " is not a PionMinus!" << G4endl;
      G4Exception("FAIL: *** exit program ***\n");
      //     return ;
    }

  // Test the DoIt for the Pion Absorption

  G4Material* apttoMaterial ;
  apttoMaterial = (*theMaterialTable)(imat) ;
  LogicalFrame->SetMaterial(apttoMaterial); 
  aPionMinus.SetKineticEnergy(0.*MeV);
  aPionMinus.SetMomentumDirection(0., 0., 1.);
  G4VParticleChange* aParticleChange;
  G4Track* aFinalParticle;
  G4String aParticleName;  
  
  thePiAbsorptionProcess->SetVerboseLevel(verboseLevel);
  
  G4int iteration = 0;   

  do {
    
    if (processId == 1)
      { aParticleChange = thePiAbsorptionProcess->AtRestDoIt(aTrack, aStep); }
    else 
      {
	if (processId == 2) 
	  { aParticleChange = gheishaAbsorptionProcess->AtRestDoIt(aTrack, aStep); }
      }
    
    // Loop over final particle List

    G4double e = 0;
    G4double eKin = 0;
    G4double Px = 0;
    G4double Py = 0;
    G4double Pz = 0;


    hNSec->accumulate(aParticleChange->GetNumberOfSecondaries());
    hDebug->accumulate(aParticleChange->GetLocalEnergyDeposit());

    G4int i;
    for (i=0; i<(aParticleChange->GetNumberOfSecondaries()); i++) 
      {
	// The following two items should be filled per event, not
	// per secondary; filled here just for convenience, to avoid
	// complicated logic 
	ntuple->column("nsec",aParticleChange->GetNumberOfSecondaries());
	ntuple->column("excitation",aParticleChange->GetLocalEnergyDeposit());

	aFinalParticle = aParticleChange->GetSecondary(i) ;
      
	e    =  aFinalParticle->GetTotalEnergy();
	eKin = aFinalParticle->GetKineticEnergy();
	Px   = (aFinalParticle->GetMomentum()).x() ;
	Py   = (aFinalParticle->GetMomentum()).y() ;
	Pz   = (aFinalParticle->GetMomentum()).z() ;
      
	aParticleName = aFinalParticle->GetDefinition()->GetParticleName();

	if (aFinalParticle->GetDefinition() == G4Proton::ProtonDefinition() ||
	    aFinalParticle->GetDefinition() == G4Neutron::NeutronDefinition())
	  {
	    hEKin->accumulate(eKin);
	    hP->accumulate(std::sqrt(Px*Px+Py*Py+Pz*Pz));

	    //	ntuple->column("px", Px);
	    //	ntuple->column("py", Py);
	    //	ntuple->column("pz", Pz);
	    //        ntuple->column("p", sqrt(Px*Px+Py*Py+Pz*Pz));
	    //        ntuple->column("e", e);
	    if (hnt == 1)
	      {
		ntuple->column("ekin", eKin);
		ntuple->dumpData(); 
	      }

	    delete aParticleChange->GetSecondary(i);
	    
	  }
      }

    G4cout << "******* Iteration = " << iteration << G4endl;
    iteration++;
    aParticleChange->Clear();
    
  }  while (iteration < nIter) ; // end of do-while
  
  hbookManager->write();

  // Clean up
  //  delete aFinalParticle;
  delete Al;
  delete Fe;
  delete Pb;
  delete Graphite;
  //  delete H;
  delete O;
  //  delete C;
  //  delete Cs;
  //  delete I;
  //  delete water;
  //  delete ethane;
  //  delete csi;

  delete touche;
  delete Step;

  
  return EXIT_SUCCESS;
}
