// -------- Start run processing

    G4StateManager* g4State=G4StateManager::GetStateManager();
    if (! g4State->SetNewState(G4State_Init)) {
      G4cout << "error changing G4state"<< G4endl;;   
    }
    G4cout << "###### Start new run # " << run << "   for "
	   << nevt << " events  #####" << G4endl;
    
    // -------- Target 

    material = mate->GetMaterial(nameMat);
    if(!material) {
      G4cout << "Material <" << nameMat
	     << "> is not found out"
	     << G4endl;
	     exit(1);
    }
    const G4Element* elm = material->GetElement(0); 

// Uzhi 1.02.13 -----------------------
    G4ProductionCuts* cuts = new G4ProductionCuts();
    G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material,cuts);
    couple->SetIndex(0);
    G4ProductionCutsTable* pkt = G4ProductionCutsTable::GetProductionCutsTable();
    std::vector<G4double>* vc = 
      const_cast<std::vector<G4double>*>(pkt->GetEnergyCutsVector(3));
    vc->push_back(0.0); 
// Uzhi 1.02.13 -----------------------

    // ------- Define target A
    G4int A = (G4int)(elm->GetN()+0.5);
    G4int Z = (G4int)(elm->GetZ()+0.5);

    // -------- Projectile
    G4ParticleDefinition* part(0);
    if (!ionParticle) {
      part = (G4ParticleTable::GetParticleTable())->FindParticle(namePart);
      if(part) energy=std::sqrt(sqr(Plab)+sqr(part->GetPDGMass()))-
                           part->GetPDGMass();
    } else {
      part = (G4ParticleTable::GetParticleTable())->GetIon(ionZ, ionA, 0.);
      if(part) energy=std::sqrt(sqr(Plab)+sqr(940.*MeV))-940.*MeV;
    }
    if (!part) {
      G4cout << " Sorry, No definition for particle" <<namePart 
	     << " found" << G4endl;
    }

    if ( ionParticle ) {       // VU
       energy*=ionA;           // VU
    }

    G4DynamicParticle dParticle(part,aDirection,energy);

    // ------- Select model
    G4VProcess* proc = phys->GetProcess(nameGen, part, material);      // Uzhi 1.02.13

    if(!proc) {
      G4cout << "For particle: " << part->GetParticleName()
	     << " generator " << nameGen << " is unavailable"<< G4endl;
	     exit(1);
    }
    G4HadronicProcess* extraproc = 0;

    G4cout<<"targetA "<<targetA<<G4endl;
    if(targetA > 0) A = targetA;
    phys->SetA(targetA);

    G4cout << "The particle:  " << part->GetParticleName() << G4endl;
    if(verbose > 0) {
      G4cout << "The material:  " << material->GetName() 
	     << " Z= " << Z << " A= " << A<< G4endl;
      G4cout << "The step:      " << theStep/mm << " mm" << G4endl;
      G4cout << "The position:  " << aPosition/mm << " mm" << G4endl;
      G4cout << "The direction: " << aDirection << G4endl;
      G4cout << "The time:      " << aTime/ns << " ns" << G4endl;
    }

    G4cout <<G4endl<< "energy = " << energy/MeV << " MeV   " 
	   << " RMS(MeV)= " << sigmae/MeV << " Plab "<<Plab/GeV<<" GeV/c"<<G4endl;

    // -------- Normalisation

    G4double cross_sec = 0.0;
    G4double cross_secel = 0.0;
    G4double cross_inel = 0.0;

    if((part == anti_proton)   || (part == anti_neutron) ||
       (part == anti_deuteron) || (part == anti_triton)  ||
       (part == anti_He3)      || (part == anti_alpha)     )
    {
     G4VComponentCrossSection* cs =0;           // Uzhi
     cs=new G4ComponentAntiNuclNuclearXS();     // Uzhi

     cs->BuildPhysicsTable(*part);
     cross_sec = cs->GetTotalElementCrossSection(part, energy, Z, (G4double) A);
     cross_inel=cs->GetInelasticElementCrossSection(part, energy, Z, (G4double)A);
     cross_secel=cs->GetElasticElementCrossSection(part, energy, Z, (G4double)A);
    }
    else
    {
     G4VCrossSectionDataSet* cs = 0;

     if(nameGen == "LElastic" || nameGen == "BertiniElastic") {
       cs = new G4HadronElasticDataSet();
     } else if (nameGen == "chargeex" ||
               nameGen == "elastic" || 
	       nameGen == "HElastic" || 
	       nameGen == "DElastic") {

       if(part == proton || part == neutron) {
	 if(xsbgg) extraproc->AddDataSet(new G4BGGNucleonElasticXS(part));
       } else if(part == pip || part == pin) {
	 if(xsbgg) extraproc->AddDataSet(new G4BGGPionElasticXS(part));
       }

     } else if(part == proton && Z > 1 && nameGen != "lepar") {
       if(xsbgg) cs = new G4BGGNucleonInelasticXS(part);
       else      cs = new G4ProtonInelasticCrossSection();

     } else if(part == neutron && Z > 1 && nameGen != "lepar") {
       if(xsbgg) cs = new G4BGGNucleonInelasticXS(part);
       else      cs = new G4NeutronInelasticCrossSection();

     } else if((part == pin || part == pip) && Z > 1 && nameGen != "lepar") {
       if(xsbgg) cs = new G4BGGPionInelasticXS(part);
       else cs = new G4PiNuclearCrossSection();

     } else if( ionParticle && !Shen ) { // Uzhi
         G4double zz = G4double(Z);
         G4double aa = G4double(A);
         if(part == deu || part == tri ||part == he3 ||part == alp)
         {
          cs = new G4TripathiLightCrossSection();
          if(cs->IsIsoApplicable(&dParticle,Z,A)) { 
	   G4cout << "Using Tripathi Light Cross section for Ions" << G4endl;
          } else {
	   delete cs;
	   cs = 0;
          }
         } else {
         if(!cs) {
	   cs = new G4TripathiCrossSection();
	   if(cs->IsIsoApplicable(&dParticle,Z,A)) {
	     G4cout << "Using Tripathi Cross section for Ions" << G4endl;
	   } else {
	     delete cs;
	     cs = 0;
	   }
          }
         }

         if(!cs) {
	   cs = new G4IonsShenCrossSection();
	   if(cs->IsElementApplicable(&dParticle,Z,material)) { 
	     G4cout << "Using Shen Cross section for Ions" << G4endl;
	   } else {
	     G4cout << "ERROR: no cross section for ion Z= " 
	   	 << zz << " A= " << aa << G4endl;
	     exit(1);
	   }
         }
     } else {
       cs = new G4HadronInelasticDataSet();
       G4cout<<"Using G4HadronInelasticDataSet()"<<G4endl;
     }

     if(extraproc) {
       extraproc->PreparePhysicsTable(*part);
       extraproc->BuildPhysicsTable(*part);
       cross_sec = extraproc->GetMicroscopicCrossSection(&dParticle,
							elm, 0);     // Uzhi 0.0);

     } else if(cs) {
       cs->BuildPhysicsTable(*part);
       cross_sec = cs->GetCrossSection(&dParticle, elm);
       G4cout<<"Try 1 cross_secel "<<cross_secel/millibarn<<G4endl;
     } else {
       cross_sec = (G4HadronCrossSections::Instance())->
       GetInelasticCrossSection(&dParticle, Z,A);         
       G4cout<<"Try 2 cross_sec "<<cross_sec/millibarn<<G4endl;
     }

// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     cs->BuildPhysicsTable(*part);                          
    } // ----------------------------------------- 

    // -------- Track
    G4Track* gTrack;
    gTrack = new G4Track(&dParticle,aTime,aPosition);
    G4TouchableHandle fpTouchable(new G4TouchableHistory());
    gTrack->SetTouchableHandle(fpTouchable);

    // -------- Step

    G4Step* step;
    step = new G4Step();
    step->SetTrack(gTrack);
    gTrack->SetStep(step);
    
    G4StepPoint *aPoint, *bPoint;
    aPoint = new G4StepPoint();
    aPoint->SetPosition(aPosition);
    aPoint->SetMaterial(material);
    G4double safety = 10000.*cm;
    aPoint->SetSafety(safety);
    step->SetPreStepPoint(aPoint);

    bPoint = aPoint;
    G4ThreeVector bPosition = aDirection*theStep;
    bPosition += aPosition;
    bPoint->SetPosition(bPosition);
    step->SetPostStepPoint(bPoint);
    step->SetStepLength(theStep);

//-------------------------------------------------------
    if(!G4StateManager::GetStateManager()->SetNewState(G4State_Idle))
      G4cout << "G4StateManager PROBLEM! " << G4endl;
    G4RotationMatrix* rot = new G4RotationMatrix();
    G4double phi0 = aDirection.phi();
    G4double theta0 = aDirection.theta();
    rot->rotateZ(-phi0);
    rot->rotateY(-theta0);

    if(verbose > 0) {
      G4cout << "Test rotation= " << (*rot)*(aDirection) << G4endl;
    }
    G4Timer* timer = new G4Timer();
    timer->Start();