// %%%%%%%%%%%%%
// gemc headers
// %%%%%%%%%%%%%
#include "DC_hitprocess.h"
#include "CLHEP/Random/RandGaussT.h"

PH_output DC_HitProcess :: ProcessHit(MHit* aHit, gemc_opts Opt)
{
 string hd_msg        = Opt.args["LOG_MSG"].args + " DC Hit Process " ;
 double HIT_VERBOSITY = Opt.args["HIT_VERBOSITY"].arg;

 PH_output out;
 out.identity = aHit->GetId();

 HCname = "DC Hit Process";
 // %%%%%%%%%%%%%%%%%%%
 // Raw hit information
 // %%%%%%%%%%%%%%%%%%%
 int nsteps = aHit->GetPos().size();



 // Identifying the fastest track in the hit
 int trackId    = -1;
 double minTime = 10000;
 vector<int>    stepTrackId = aHit->GetTrackId();
 vector<double> stepTime    = aHit->GetTime();
 vector<G4double> Edep = aHit->GetEdep();

 for(int s=0; s<nsteps; s++)
 {
    if(stepTrackId[s] != trackId && stepTime[s]/ns < minTime && Edep[s] >= aHit->GetThreshold())
       trackId = stepTrackId[s];
 }

 // If no step pass the threshold, getting the fastest track
 if(trackId == -1)
 {
    for(int s=0; s<nsteps; s++)
    {
       if(stepTrackId[s] != trackId && stepTime[s]/ns < minTime)
          trackId = stepTrackId[s];
    }
 }

 // Get Total Energy deposited of the fastest track
 double Etot = 0;
 int count_track_steps = 0;
 for(int s=0; s<nsteps; s++)
    if(stepTrackId[s] == trackId)
    {
       Etot = Etot + Edep[s];
       count_track_steps++;
    }


 // average global, local positions of fastest track
 double x, y, z;
 double lx, ly, lz;
 x = y = z = lx = ly = lz = 0;
 vector<G4ThreeVector> pos  = aHit->GetPos();
 vector<G4ThreeVector> Lpos = aHit->GetLPos();

/* if(Etot>0)*/
 for(int s=0; s<nsteps; s++)
 {
    if(stepTrackId[s] == trackId)
    {
/*       x  = x  +  pos[s].x()*Edep[s]/Etot;
       y  = y  +  pos[s].y()*Edep[s]/Etot;
       z  = z  +  pos[s].z()*Edep[s]/Etot;
       lx = lx + Lpos[s].x()*Edep[s]/Etot;
       ly = ly + Lpos[s].y()*Edep[s]/Etot;
       lz = lz + Lpos[s].z()*Edep[s]/Etot;*/
       x  = x  +  pos[s].x()/count_track_steps;
       y  = y  +  pos[s].y()/count_track_steps;
       z  = z  +  pos[s].z()/count_track_steps;
       lx = lx + Lpos[s].x()/count_track_steps;
       ly = ly + Lpos[s].y()/count_track_steps;
       lz = lz + Lpos[s].z()/count_track_steps;
    }
 }
//  else
//  {
//    x  = pos[0].x();
//    y  = pos[0].y();
//    z  = pos[0].z();
//    lx = Lpos[0].x();
//    ly = Lpos[0].y();
//    lz = Lpos[0].z();
//  }

 // average time of fastest track
 double time = 0;
 vector<G4double> times = aHit->GetTime();
 for(int s=0; s<nsteps; s++) if(stepTrackId[s] == trackId) time = time + times[s]/nsteps;

 // Energy of the track
 double Ene = aHit->GetE();

 out.raws.push_back(Etot);
 out.raws.push_back(x);
 out.raws.push_back(y);
 out.raws.push_back(z);
 out.raws.push_back(lx);
 out.raws.push_back(ly);
 out.raws.push_back(lz);
 out.raws.push_back(time);
 out.raws.push_back((double) aHit->GetPID());
 out.raws.push_back(aHit->GetVert().getX());
 out.raws.push_back(aHit->GetVert().getY());
 out.raws.push_back(aHit->GetVert().getZ());
 out.raws.push_back(Ene);
 out.raws.push_back((double) aHit->GetmPID());
 out.raws.push_back(aHit->GetmVert().getX());
 out.raws.push_back(aHit->GetmVert().getY());
 out.raws.push_back(aHit->GetmVert().getZ());



 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Distance 1: smear by 300 um
 // Drift Time 1: using 50 um/ns
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%

 // Finding Wire position
 double NWIRES = 112;
 double ylength =  aHit->GetDetector().dimensions[3];  
 double deltay  = 2.0*ylength/NWIRES;
 vector<identifier> iden = aHit->GetId();
 int nwire      = iden[3].id;
 double WIRE_Y  = nwire*deltay;

 // Finding DOCA
 double doca    = 10000;
 double LR      = 0;

 for(int s=0; s<nsteps; s++)
 {
   G4ThreeVector DOCA(0, Lpos[s].y() + ylength - WIRE_Y, Lpos[s].z());
   if(DOCA.mag() <= doca && stepTrackId[s] == trackId )
   {
      doca = DOCA.mag();
      if(DOCA.y() >=0 ) LR = 1;
      else  LR = -1;

   }
 }
//   if(aHit->GetPID() == 2112)
//   cout <<  "    wire=" << nwire << "    layer=" << iden[2].id  
//        <<  "    doca=" << doca << "     nsteps=" << count_track_steps << endl;




 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Digitization
 // DC ID:
 // sector, SL, Layer, Wire
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 int sector = out.identity[0].id;
 int SL     = out.identity[1].id;
 int Layer  = out.identity[2].id;
 int Wire   = out.identity[3].id;

 double drift_velocity = 0;
 if(SL == 1 || SL == 2) drift_velocity = 0.053;    
 if(SL == 3 || SL == 4) drift_velocity = 0.026;    
 if(SL == 5 || SL == 6) drift_velocity = 0.036;    

 double sdoca  = -1;
 while(sdoca < 0) sdoca = CLHEP::RandGauss::shoot(doca, 0.3);  
 double time1  = doca/drift_velocity;
 double stime1 = sdoca/drift_velocity;

 out.raws.push_back(LR);
 out.raws.push_back(doca);
 out.raws.push_back(sdoca);
 out.raws.push_back(time1);
 out.raws.push_back(stime1);


 out.dgtz.push_back(sector);
 out.dgtz.push_back(SL);
 out.dgtz.push_back(Layer);
 out.dgtz.push_back(Wire);

 if(HIT_VERBOSITY>4)
 cout <<  hd_msg << "  Sector: "     << sector
                 << "  Superlayer: " << SL
                 << "  Layer: "      << Layer
                 << "  Wire: "       << Wire
                 << "  L/R: "        << LR
                 << "  doca: "       << doca
                 << "  sdoca: "      << sdoca
                 << "  time1"        << time1
                 << "  stime1="      << stime1 <<  endl;

 return out;
}

vector<identifier>  DC_HitProcess :: ProcessID(vector<identifier> id, G4Step* aStep, detector Detector)
{
 double NWIRES = 112;

 G4StepPoint   *prestep   = aStep->GetPreStepPoint();
 G4StepPoint   *poststep  = aStep->GetPostStepPoint();
 G4VTouchable* TH = (G4VTouchable*) aStep->GetPreStepPoint()->GetTouchable();

 string         name    = TH->GetVolume(0)->GetName();                                    
 G4ThreeVector   xyz    = poststep->GetPosition();                                        
 G4ThreeVector  Lxyz    = prestep->GetTouchableHandle()->GetHistory()                     
                         ->GetTopTransform().TransformPoint(xyz);

 double ylength = Detector.dimensions[3];  
 double deltay  = 2.0*ylength/NWIRES;
 double loc_y   = Lxyz.y() + ylength;      

 int nwire = (int) floor(loc_y/deltay);
 vector<identifier> yid = id;

 if(nwire > NWIRES) cout << " SuperLayer: " << yid[1].id << "      layer: "        << yid[2].id
                         << "       wire: " << nwire     << "      length: "       << ylength
                         << "         ly: " << Lxyz.y()  << "      deltay: "       << deltay
                         << "      loc_y: " << loc_y     << "      nwire*deltay: " << fabs( nwire*deltay - loc_y ) << endl;


 yid[3].id = nwire;

 if(fabs( (nwire+1)*deltay - loc_y ) < fabs( nwire*deltay - loc_y )  )
    yid[3].id = nwire + 1;

 return yid;
}