// $Id: PamVMCPrimaryGenerator.h,v 1.0 2007/06/03 

#ifndef PAMVMC_PRIMARY_GENERATOR_H
#define PAMVMC_PRIMARY_GENERATOR_H

#include <iostream>
#include <TVirtualMCApplication.h>
#include <TClonesArray.h>
#include <TVector3.h>
#include <TVector2.h>
#include <TMath.h>
#include <TRandom.h>
#include <TF1.h>
//#include <TObjectTable.h>

#include "PamRootManager.h"
#include "PamVMC2PiAux.h"
#include "PamVMCPrimaryGF.h"

class TVirtualMCStack;

using std::cout;
using std::endl;
using TMath::Sqrt;

class PamVMCPrimary : public TObject
{

 public:
  
 PamVMCPrimary(): fID(0), fPDG(0), fX0(0.), fY0(0.), fZ0(0.), fTHETA(0.), fPHI(0.), fP0(0.), fGOOD(kFALSE),  fXTOL(-100.), fYTOL(-100.){
    

    for (Int_t j = 0; j<14; j++){ 
      fCR_X[j] = -100.;	
      fCR_Y[j] = -100.;		
    }	    
  };	

   virtual ~PamVMCPrimary(){ /*gObjectTable->Remove(this);*/};

  UInt_t fID;
  Int_t fPDG;
  Double_t fX0, fY0, fZ0;
  Double_t fTHETA, fPHI;
  Double_t fP0;
  Bool_t fGOOD;
  Double_t fXTOL, fYTOL;
  Double_t fCR_X [14];
  Double_t fCR_Y [14];
  

  void Clean() { fID=0; fPDG=0; fX0=fY0=fZ0=fTHETA=fPHI=fP0=0.; fGOOD=kFALSE; fXTOL=fYTOL=-100.;
    for (Int_t j = 0; j<14; j++){
      fCR_X[j] = -100.;
      fCR_Y[j] = -100.;      
    }	
    
  }

  void Clear(const Option_t * =""){ this->~PamVMCPrimary(); }

  void Print(const Option_t * ="") const
    { cout<<"PRIMARY particle "<<fID<<"  infromation:"<<endl;
      cout<<"Pdg="<<fPDG<<endl; 
      cout<<"Position: "<<"("<<fX0<<","<<fY0<<","<<fZ0<<")"<<endl;
      cout<<"P0, Theta, Phi: "<<fP0<<","<<fTHETA<<","<<fPHI<<endl;
      cout<<"GOOD Single Track:"<<fGOOD<<endl;
      cout<<"Xtol= "<<fXTOL<<" Ytol= "<<fYTOL<<endl;
      cout<<" Coordinates of acceptance planes crossing:"<<endl;
      for(Int_t i = 0; i< 14; i++){
	cout<<"i="<<i<<" X="<<fCR_X[i]<<" Y="<<fCR_Y[i]<<endl;
      }
    }

  ClassDef(PamVMCPrimary,1);
};

PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b);

class PamVMCPrimaryGenerator : public TObject
{
  public:
    PamVMCPrimaryGenerator(TVirtualMCStack* stack); 
    PamVMCPrimaryGenerator();

    virtual ~PamVMCPrimaryGenerator();

    // methods
    
    void GeneratePrimary();
    // set methods
    void SetParticle(Int_t pdg);
    void SetParticleID(UInt_t id){ fprim.fID = id; };
    void SetPosition(Double_t vx, Double_t vy, Double_t vz) 
         {fprim.fX0=vx; fprim.fY0=vy; fprim.fZ0=vz;};
   
    void SetKinEnergy(Double_t kinEnergy) {fprim.fP0=KinEToMomentum(kinEnergy);};
    void SetRigidity(Double_t rigidity) {fprim.fP0=RigToMomentum(rigidity);};
    void SetMomentum(Double_t momentum) {fprim.fP0=momentum; };

    void SetDirection(Double_t theta, Double_t phi) {fprim.fTHETA=theta; fprim.fPHI=phi; };
    void SetMomentum(Double_t px, Double_t py, Double_t pz);
    
     
    // gen methods

    void GenPosition(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax){
      SetPosition(frandom->Uniform(xmin,xmax),frandom->Uniform(ymin,ymax),frandom->Uniform(zmin,zmax)); }
    
    void GenDirection(Double_t thetamin, Double_t thetamax, Double_t phimin, Double_t phimax){
      Double_t theta = TMath::ASin(Sqrt( frandom->Uniform( TMath::ASin(thetamin), TMath::ASin(thetamax))));
      Double_t phi = TMath::TwoPi() * frandom->Uniform( phimin, phimax );
      SetDirection(theta, phi);
    }
    
    void GenPosAng(TVector3 &pos, TVector2 &ang) { f2PiAux->GenPosAng(pos,ang); }
    
    void Gen2PiPosAng(){
      TVector3 pos;
      TVector2 ang;
      GenPosAng(pos,ang);
      SetPosition( pos.X(), pos.Y(), pos.Z() );
      SetDirection( ang.X(), ang.Y() );
    }
    
    //flat spectra generator
    void GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy=kFALSE); 
    //power law spectra, gamma - differential spectral index
    void GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE); 
    void GenSpe_Flat(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE); 
    void GenSpe_3par(Double_t PEmin, Double_t PEmax, Double_t a, Double_t b, Double_t c);
    void GenSphericalPhiThe(); // cecilia // flusso istropo // sets phi the
    void GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, 
		      Double_t zmin, Double_t zmax); // flusso isotropo entro S2 S3
    //                                                  sets position and phi the
    

    // get methods
    Int_t GetParticleID(){ return fprim.fID; };
    Int_t GetParticle(){ return fprim.fPDG; };
    void  GetPositon(TVector3 & v){ v.SetXYZ(fprim.fX0,fprim.fY0,fprim.fZ0); };
    void  GetPosition(Double_t & X0, Double_t & Y0,Double_t & Z0)
          {X0=fprim.fX0; Y0=fprim.fY0; Z0=fprim.fZ0;};

    Double_t GetKinEnergy() { return MomentumToKinE(fprim.fP0); };
    Double_t GetRigidity() { return MomentumToRig(fprim.fP0); };
    Double_t GetMomentum() { return fprim.fP0; };
    
    Double_t Get2PiArea() { return f2PiAux->GetArea(); };

    void GetDirection(TVector2 & v) {v.Set(fprim.fTHETA,fprim.fPHI);}
    void GetDirection(Double_t &theta, Double_t & phi) { theta = fprim.fTHETA; phi = fprim.fPHI;}
    
    Bool_t Getgood(){ return fprim.fGOOD; };

    void GetTol(Double_t& xtol, Double_t& ytol){ xtol=fprim.fXTOL; ytol=fprim.fYTOL; };

    void GetCR_GF(Double_t& xcr, Double_t& ycr, Int_t id_GF){
      if(id_GF<14){
	xcr=fprim.fCR_X[id_GF]; 
	ycr=fprim.fCR_Y[id_GF];
      }
    }
    
    Int_t GetNgood(){ return fngood; }

    //initialize random
    void SetRandom(TRandom* random){
      frandom = random;
      f2PiAux->SetRandom(random);
    }

    //work with collection of primaries

    void Register(){
      PamRootManager::Instance()->
	Register("PRIM","TClonesArray", &fprimColl);
    }

    void ClearPrimCol(){
      fevno=0;
      (*fprimColl).Clear("C");
    }
    //all methods to operate with PamVMCPrimaryGF object
    //after each step of primary
    void SetZeroStep(Double_t z_curr, Double_t x_curr, Double_t y_curr){
      fzprev = z_curr;
      fxprev = x_curr;
      fyprev = y_curr;
    }
    void CheckCrossing(Double_t z_curr, Double_t x_curr, Double_t y_curr){
      fpgf->CheckCrossing(fzprev,fxprev,fyprev, z_curr, x_curr, y_curr);
      SetZeroStep(z_curr,x_curr, y_curr);
    }
    //after last step of primary
    void CheckInsideAcceptance(){
      for(Int_t i = 0; i<NGF; i++) fprim.fCR_X[i] = fprim.fCR_Y[i] = fprim.fXTOL = fprim.fYTOL = -100.; //zeroing
      fprim.fGOOD = fpgf->IsInsideAcceptance(fprim.fCR_X, fprim.fCR_Y, fprim.fXTOL, fprim.fYTOL);
      fpgf->Clear();
      if(fprim.fGOOD) fngood++;
      //Here I'm getting what I have in collection and asume that I have only 1 primary here
      //If I wish to launch more in 1 event I need change this stuff
      PamVMCPrimary* tmp = (PamVMCPrimary*)fprimColl->At(0);
      if(tmp) *tmp = fprim;
    }

    Double_t MomentumToKinE(Double_t P0);
    Double_t KinEToMomentum(Double_t E0);
    Double_t MomentumToRig(Double_t P0) { return TMath::Abs(P0/fcharge); };
    Double_t RigToMomentum(Double_t R0){ return TMath::Abs(R0*fcharge); };
    Double_t function3par(Double_t xx, Double_t a, Double_t b, Double_t c){return  5.*pow((xx + b * exp(-c * sqrt(xx))),-a);};

  private:
    // methods
    
  

    // data members
    TVirtualMCStack*  fStack;
    Int_t             fevno;
    Int_t             fngood; //particles inside nominal acceptance
    PamVMCPrimary     fprim;
    Double_t          fmass;
    Double_t          fcharge;
    Double_t          fxprev, fyprev, fzprev; // coordinates of primary on prev. step
    TClonesArray*     fprimColl; 
    TRandom*          frandom; // Class is not a owner of this object
    TF1*              ftheta; // To generate sherical distributhin in theta-angle
    PamVMC2PiAux*     f2PiAux; //2 Pi generator aux stuff
    PamVMCPrimaryGF*  fpgf;  // Object which
    

  ClassDef(PamVMCPrimaryGenerator,1)  //PamVMCPrimaryGenerator
};

// inline functions


inline Double_t PamVMCPrimaryGenerator::MomentumToKinE(Double_t P0)
{ return  Sqrt(P0*P0+fmass*fmass)-fmass; }

inline Double_t PamVMCPrimaryGenerator::KinEToMomentum(Double_t E0)
{ return  Sqrt(E0*E0+2.*E0*fmass); }

#endif //PAMVMC_PRIMARY_GENERATOR_H