#include <vector>

#include <RAT/GenPMTEff.hh>
#include <RAT/Log.hh>
#include <RAT/GLG4StringUtil.hh>
#include <RAT/GLG4TimeGen.hh>
#include <RAT/DB.hh>
#include <RAT/DU/ChanHWStatus.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/Utility.hh>
#include <Randomize.hh>

#include <G4Event.hh>
#include <G4OpticalPhoton.hh>
#include <G4PhysicalConstants.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4ThreeVector.hh>


using CLHEP::nanometer;
using CLHEP::twopi;
using namespace std;

namespace RAT {


GenPMTEff::GenPMTEff() : GLG4Gen("PMTEff"), fStateStr("") {
    fTimeGen = new GLG4TimeGen_Poisson();
    fPrintNOnLine = true;
}

GenPMTEff::~GenPMTEff() {
}


void GenPMTEff::BeginOfRun() {
    // Initialize underlying low level generators
    GLG4Gen::BeginOfRun();
    DB* db = DB::Get();

    DBLinkPtr ptrInnerPMT = db->GetLink( "GEO", "innerPMT" );
    vector<string> conType = ptrInnerPMT->GetSArray("concentrator_type");
    DBLinkPtr ptrConcentrator = db->GetLink( "CONCENTRATOR", conType[0]);
    rConc = ptrConcentrator->GetDArray("rho_coord").back() -0.1;
}

void GenPMTEff::GenerateEvent(G4Event* event) {
    G4double t0 = NextTime();
    G4ThreeVector pos;

    DU::PMTInfo pmtInfo = DU::Utility::Get()->GetPMTInfo();
    const DU::ChanHWStatus&
                channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();
    size_t nPMTs = pmtInfo.GetCount();


    G4PrimaryParticle* primary;
    int nOnLine = 0;
    for(int k=0; k<(int)nPMTs; k++) {
        bool isOnLine = (pmtInfo.GetType(k))== pmtInfo.NORMAL ;
        if( channelHardwareStatus.IsEnabled() ) {       //  Skip off-line PMTs
            isOnLine = isOnLine && channelHardwareStatus.IsTubeOnline(k)
                                && channelHardwareStatus.IsChannelOnline(k);
        }
        if(isOnLine) {
            nOnLine++;
            // Vector of PMT location (actually is center of concentrator
            // "face"):
            TVector3 posPMT = pmtInfo.GetPosition(k);
            // dirPMT is PMT orientation vector, the direction from
            // concentrator center towards the PMT itself
            TVector3 dirPMT = pmtInfo.GetDirection(k);

            // Arbitrary unit vector perpendicular to PMT orientation vector:
            TVector3 ranPerp = (dirPMT.Orthogonal()).Unit();
            // Rotate ranPerp by random angle around PMT orientation vector to
            // get random unit vector perpendicular to PMT orientation vector:
            ranPerp.Rotate( twopi*G4UniformRand(), dirPMT );
            // Draw random distance away from PMT axis, up to maximum of rConc,
            // such that one will have a uniform distribution in area over the
            // concentrator face:
            double rho = rConc*sqrt(G4UniformRand());
            double phi = twopi*G4UniformRand();
            // position photon uniformly on concentrator face:
            TVector3 posPhoton = posPMT + cos(phi)*rho*ranPerp
                                        + sin(phi)*rho*dirPMT.Cross(ranPerp) ;
            // Another random axis perpendicular to PMT orientation vector:
            ranPerp.Rotate( twopi*G4UniformRand(), dirPMT );
            // Set photon direction initially equal to PMT orientation, i.e.,
            // outward perpendicular to the concentrator face:
            TVector3 photonMom = fPhEnergy*dirPMT;
            // Now rotate about randomly chosen axis by desired incidence angle:
            photonMom.Rotate( fIncAngle, ranPerp );
            G4ThreeVector g4PhotonMom =
                    G4ThreeVector(photonMom.X(), photonMom.Y(), photonMom.Z());
            primary = new G4PrimaryParticle( G4OpticalPhoton::Definition(),
                  g4PhotonMom.x(), g4PhotonMom.y(), g4PhotonMom.z());

            G4ThreeVector pol;
            if ( fIncAngle==0.) {
                // Generate random polarization:
                phi = twopi* G4UniformRand();
                G4ThreeVector e1 = g4PhotonMom.orthogonal().unit();
                G4ThreeVector e2 = g4PhotonMom.unit().cross(e1);
                pol = e1*cos(phi)+e2*sin(phi);
            }else{
                // Generate desired polarization
                // Unit vector perpendicular to plane of incidence:
                TVector3 unitSenkrecht = (dirPMT.Cross(photonMom)).Unit();
                if( fPSenkrecht ) {     // polarization perpendicular to plane
                    pol = G4ThreeVector( unitSenkrecht.X(), unitSenkrecht.Y(),
                                         unitSenkrecht.Z()) ;
                }else{                  // polarization in the plane
                    TVector3 unitParallel =
                                        (unitSenkrecht.Cross(photonMom)).Unit();
                    pol = G4ThreeVector( unitParallel.X(), unitParallel.Y(),
                                         unitParallel.Z()) ;
                }
            }
            primary->SetPolarization(pol);

            G4ThreeVector g4PhotonVtx =
                    G4ThreeVector(posPhoton.X(), posPhoton.Y(), posPhoton.Z());
            G4PrimaryVertex* vertex = new G4PrimaryVertex(g4PhotonVtx, t0);
            vertex->SetPrimary(primary);
            event->AddPrimaryVertex(vertex);
        }       // if(isOnLine)
    }           // for(int k=0; k<(int)nPMTs; k++)
    if(fPrintNOnLine) {
        fPrintNOnLine = false;
        info << "GenPMTEff:  Number of tubes on-line = " << nOnLine << endl;
    }
}

void GenPMTEff::ResetTime(double offset) {
    fNextTime = fTimeGen->GenerateEventTime() + offset;
}

void GenPMTEff::SetState(G4String state) {
    fStateStr = state;
    // arguments are WL:IA:PP, where WL=wavelength, IA=incidence angle, and
    // PP=polarization
    // lambda = desired photon wavelength in nm
    // fIncAngle = incidence angle of photons on concentrator face, in degrees
    // fPSenkrect = true for senkrecht, false for parallel
    state = util_strip_default(state);
    vector<string> parts = util_split(state, ":");
    // Defaults:
    double lambda = 385.;
    fIncAngle = 0.;
    fPSenkrecht = true;

    for(unsigned int i=0; i<parts.size(); i++) {
        if(i==0) lambda = util_to_double(parts[i]);
        if(i==1) fIncAngle = util_to_double(parts[i]);
        if(i==2) fPSenkrecht = ( (parts[i])[0]=='s' );
    }
    fPhEnergy = twopi*hbarc/(lambda*nanometer);
    info << "GenPMTEff: lambda, fPhEnergy, fIncAngle, fPSenkrecht = "
        << lambda <<" "<< fPhEnergy <<" "<< fIncAngle <<" "<<
        fPSenkrecht << endl;
    fIncAngle = fIncAngle*pi/180.;      // convert angle to radians
}

G4String GenPMTEff::GetState() const {
    return fStateStr;
}

} // namespace RAT