//
// C++ Implementation: PndMvdCalcStrip
//
// Description: 
//
//
// Author: HG Zaunick <hg.zaunick@physik.tu-dresden.de>, (C) 2007
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include <cmath>

#include "PndMvdCalcStrip.h"
#include "TRandom3.h"

PndMvdCalcStrip::PndMvdCalcStrip(){
  fPitch = 0.;
  fOrient = 0.;
  fAnchor = TVector2(0.,0.);
  fNrStrips = 0;
  fThreshold = 0.;
  fNoise = 0.;
  fVerboseLevel = 1;
  fRNG = new TRandom3();
}

PndMvdCalcStrip::PndMvdCalcStrip(Double_t pitch, Double_t orient,
             Int_t nrStrips, Int_t nrFeChannels,
             const TVector2& firstStripAnchor,
             Double_t threshold, Double_t noise)
   : fPitch(pitch), fOrient(orient),
     fNrStrips(nrStrips), fNrFeChannels(nrFeChannels),
     fAnchor(firstStripAnchor),
     fThreshold(threshold), fNoise(noise)
{
      fStripDir.Set(cos(fOrient),sin(fOrient));
      fOrthoDir.Set(sin(fOrient),-cos(fOrient));
      fVerboseLevel = 1;
  fRNG = new TRandom3();
  //Print();
}



std::vector<PndMvdStrip>
PndMvdCalcStrip::GetStrips(Double_t inx, Double_t iny, Double_t inz,
                        Double_t outx, Double_t outy, Double_t outz,
                        Double_t eLoss)
{
    if (fVerboseLevel > 2) std::cout<<"-I-  PndMvdCalcStrip::GetStrips "<<std::endl;
    
    // 2d-Projection of trajectory
    TVector2 in(inx,iny);
    TVector2 out(outx,outy);
    TVector2 path=out-in;
    
    if (fVerboseLevel > 2){
        std::cout<<" InPoint: ("<<in.X()<<","<<in.Y()<<std::endl;
        std::cout<<" OutPoint: ("<<out.X()<<","<<out.Y()<<std::endl;
    }
    
    std::vector<PndMvdStrip> strips;
    Double_t fSmeared;
        
    if (path.Mod()<1E-18) {
        std::cout<<"-W-  PndMvdCalcStrip::GetStrips : No Trajectory inside Sensor!"<<path.Mod()<<std::endl;
        return strips;
    }    
    
    if (fVerboseLevel > 1) std::cout<<" pathlength: "<<path.Mod()<<std::endl;
    
    Double_t nuIn  = CalcStripFromPoint(inx,iny);
    Double_t nuOut = CalcStripFromPoint(outx,outy);
    
    if (fVerboseLevel > 2) std::cout<<" nuIn = "<<nuIn<<" ; nuOut = "<<nuOut<<std::endl;
    
    Double_t dir = 0.;
    if (nuOut>nuIn) dir=1.;
    else dir = -1.;
    
    Double_t Q = eLoss/3.61; // 3.6 eV/Electron in Silicon
    if (fVerboseLevel > 1) std::cout<<" integral charge = "<<Q<<std::endl;
    
    // did we hit the active area ?
    if ( (nuIn<0.5 && nuOut<0.5) || (((nuIn+0.5) > Double_t(fNrStrips-1)) && ((nuOut+0.5) > Double_t(fNrStrips-1)))){
        if (fVerboseLevel > 1) std::cout<<" Hit outside active area."<<std::endl;
        return strips;
    }
    
    // is the In-Point inside active area ?
    if (nuIn<-0.5){
        Q *= (nuOut+0.5)/(nuOut-nuIn);
        nuIn = -0.5;
    } else if (nuIn > (Double_t(fNrStrips)-0.5)){
        Q *= ((Double_t)fNrStrips-0.5-nuOut)/(-nuOut+nuIn);
        nuIn = Double_t(fNrStrips)-.5;
    }
    
    // is the Out-Point inside active area ?
    if (nuOut<-0.5){
        Q *= (nuIn+0.5)/(-nuOut+nuIn);
        nuOut = -0.5;
    } else if (nuOut > (Double_t(fNrStrips)-0.5)){
        Q *= ((Double_t)fNrStrips-0.5-nuIn)/(nuOut-nuIn);
        nuOut = Double_t(fNrStrips)-.5;
    }
    
    // only one strip hit ?
    if (Int_t(nuIn+0.5) == Int_t(nuOut+0.5)){
        // this strip collected the entire charge
        fSmeared = SmearCharge(Q);
        if (fSmeared >= fThreshold)
            strips.push_back(PndMvdStrip(Int_t(nuOut+0.5),fSmeared));
        if (fVerboseLevel > 1) std::cout<<" -> 1 strip hit."<<std::endl;
        return strips;
    } else {
        // calculate portion of track in first strip
        Int_t nrHits = 0;
        if (dir>0.){
            Double_t Q1 = Q*(std::floor(nuIn+0.5)+0.5-nuIn)/(nuOut-nuIn);
            if (fVerboseLevel > 2){
                std::cout<<" part of first strip : "<<std::floor(nuIn+0.5)+0.5-nuIn<<std::endl ;
                std::cout<<" charge : "<<Q1<<std::endl ;
                std::cout<<" next strip : "<<std::floor(nuIn+0.5)+0.5<<std::endl ;
            }              
            fSmeared = SmearCharge(Q1);
            if (fSmeared >= fThreshold)
                strips.push_back(PndMvdStrip(Int_t(nuIn+0.5),fSmeared));
            nrHits++;
            Q -= Q1;
            nuIn = std::floor(nuIn+0.5) + 0.5;
        } else {
            Double_t Q1 = Q*(nuIn-(std::floor(nuIn+0.5)-0.5))/(-nuOut+nuIn);
            if (fVerboseLevel > 2){
                std::cout<<" part of first strip : "<<nuIn-(std::floor(nuIn+0.5)-0.5)<<std::endl ;
                std::cout<<" charge : "<<Q1<<std::endl ;
                std::cout<<" next strip : "<<std::floor(nuIn+0.5)-0.5<<std::endl ;
            }
            fSmeared = SmearCharge(Q1);
            if (fSmeared >= fThreshold)
                strips.push_back(PndMvdStrip(Int_t(nuIn+0.5),fSmeared));
            nrHits++;
            Q -= Q1;
            nuIn = std::floor(nuIn+0.5) - 0.5;
        }
        
        // calculate portion of track in last strip
        if (dir<0.){
            Double_t Q2 = Q*(std::floor(nuOut+0.5)+0.5-nuOut)/(nuIn-nuOut);
            if (fVerboseLevel > 2){
                std::cout<<" part of last strip : "<<std::floor(nuOut+0.5)+0.5-nuOut<<std::endl ;
                std::cout<<" charge : "<<Q2<<std::endl ;
                std::cout<<" next strip : "<<std::floor(nuOut+0.5)+0.5<<std::endl ;
            }  
            fSmeared = SmearCharge(Q2);
            if (fSmeared >= fThreshold)
                strips.push_back(PndMvdStrip(Int_t(nuOut+0.5),fSmeared));
            nrHits++;
            Q -= Q2;
            nuOut = std::floor(nuOut+0.5) + 0.5;
        } else {
            Double_t Q2 = Q*(nuOut-(std::floor(nuOut+0.5)-0.5))/(-nuIn+nuOut);
            if (fVerboseLevel > 2){
                std::cout<<" part of last strip : "<<nuOut-(std::floor(nuOut+0.5)-0.5)<<std::endl ;
                std::cout<<" charge : "<<Q2<<std::endl ;
                std::cout<<" next strip : "<<std::floor(nuOut+0.5)-0.5<<std::endl ;
            }   
            fSmeared = SmearCharge(Q2);
            if (fSmeared >= fThreshold)
                strips.push_back(PndMvdStrip(Int_t(nuOut+0.5),fSmeared));
            nrHits++;
            Q -= Q2;
            nuOut = std::floor(nuOut+0.5) - 0.5;
        }
        
        if (fVerboseLevel > 2) {
            std::cout<<" dir="<<Int_t((dir>0.)?1:-1)<<std::endl;
            std::cout<<" begin="<<Int_t(nuIn)<<" end="<<Int_t(nuOut)<<std::endl;
        }

        // Distribute the charge among the intermediate strips
        if (Int_t(nuOut)!=Int_t(nuIn)) {
            Double_t dQ=Q/std::fabs(nuOut-nuIn);
            for (Int_t n=Int_t(nuIn)+Int_t((dir>0.)?1:0); n!=Int_t(nuOut)+Int_t((dir>0.)?1:0); n+=Int_t((dir>0.)?1:-1))
            {
                if (fVerboseLevel > 3) std::cout<<" n = "<<n<<std::endl;
                fSmeared = SmearCharge(dQ);
                if (fSmeared >= fThreshold)
                    strips.push_back(PndMvdStrip(n,fSmeared));
                nrHits++;
                Q -= dQ;
            }
        }
        if (fVerboseLevel > 2) if (fabs(Q)>1.) std::cout<<" charge Q = "<<Q<<" not detected!"<<std::endl;
        if (fVerboseLevel > 1) std::cout<<" -> "<<nrHits<<" strips hit."<<std::endl;
    }
    
    return strips;
}

Double_t PndMvdCalcStrip::SmearCharge(Double_t charge)
{
    Double_t smeared = fRNG->Gaus(charge,fNoise);
    if (fVerboseLevel > 3) std::cout<<" charge = "<<charge<<", smeared = "<<smeared<<std::endl;
    return smeared;
}

Double_t PndMvdCalcStrip::CalcStripFromPoint(Double_t x, Double_t y)
{
    return ( (x-fAnchor.X())*fStripDir.Y()
             -(y-fAnchor.Y())*fStripDir.X() )/fPitch;
}

Int_t PndMvdCalcStrip::CalcFEfromStrip(Int_t stripNr) const { return (stripNr/fNrFeChannels); }

Int_t PndMvdCalcStrip::CalcChannelfromStrip(Int_t stripNr) const { return (stripNr%fNrFeChannels); }

void PndMvdCalcStrip::CalcFeChToStrip(Int_t fe, Int_t channel, Int_t& strip, enum SensorSide& side) const
{
    Int_t nr = fe * fNrFeChannels + channel;
    if (nr > fNrStrips) {
        nr -= fNrStrips;
        side = kBOTTOM;
    } else side = kTOP;
    strip = nr;
}


void PndMvdCalcStrip::CalcStripPointOnLine(Double_t strip,TVector2& point) const
{
  point =  fAnchor + ( fPitch * strip * fOrthoDir ) ;
}


void PndMvdCalcStrip::Print() const
{
    std::cout<<"-I- PndMvdCalcStrip Info :"<<std::endl;
    std::cout<<"     pitch                  = "<<fPitch*10000.<<" um"<<std::endl;
    std::cout<<"     orientation angle      = "<<fOrient/TMath::Pi()*180.<<" deg"<<std::endl;
    std::cout<<"     nr of strips           = "<<fNrStrips<<std::endl;
    std::cout<<"     nr of channels per FE  = "<<fNrFeChannels<<std::endl;
    std::cout<<"     nr of frontends        = "<<fNrStrips/fNrFeChannels<<std::endl;
    std::cout<<"     anchor point           = "<<fAnchor.X()<<","<<fAnchor.Y()<<")"<<std::endl;
    std::cout<<"     strip-direction vector = ("<<fStripDir.X()<<","<<fStripDir.Y()<<")"<<std::endl;
}

// const Double_t PndMvdCalcStrip::CalcXOnLine(const Double_t strip)
// {
//   return ( fAnchor.X() + fPitch*strip*fOrthoDir.X() );
// }
// const Double_t PndMvdCalcStrip::CalcYOnLine(const Double_t strip)
// {
//   return ( fAnchor.Y() + fPitch*strip*fOrthoDir.Y() );
// }
// 
// TVector2 PndMvdCalcStrip::CalcCrossingPoint(Double_t strip,
//            const PndMvdCalcStrip& otherside, Double_t otherstrip)
// {
// 
//   Double_t xhere = CalcXOnLine(strip);
//   Double_t yhere = CalcYOnLine(strip);
//   Double_t xthere = otherside.CalcXOnLine(strip);
//   Double_t ythere = otherside.CalcYOnLine(strip);
// 
//   Double_t dx = xthere - xhere;
//   Double_t dy = ythere - yhere;
// 
//   Double_t M = fStripDir.X()*otherside.fStripDir.Y() - fStripDir.Y()*otherside.fStripDir.X();
//   if(M != 0.)
//   {
//     Double_t s  = fStripDir.Y()*dx/M - fStripDir.X()*dy/M;
// 
//     TVector2 result(
//       otherside.fStripDir.X()*s + xthere ,
//       otherside.fStripDir.Y()*s + ythere );
//     return result;
//   }
// 
//   TVector2 r(0.,0.);
//   return r;
//   
// }