// $Header: SmearTrack.cxx $
// algorithm from K.Riles 8.12.98
// Version 1.0 12-Aug-98 Richard incorporate K.Riles' algorithm.
// Version 1.1 06-Mar-99 Richard get start point from parent MC.

#include "TString.h"
#include "SmearTrack.h"
#include "GetTrackLookups.h"
#include <math.h>
#include <stdlib.h>

//______________________________________________________________________
//
// SmearTrack
//
// Smears tracks via lookup tables on resolution      

ClassImp(SmearTrack)
 Track* SmearTrack::makeTrack(McPart* p, Int_t index) {

  // creates track for input McPart if it passes acceptance cuts and
  // uses lookup tables for resolutions

// if MCParticle pointer is zero, it belongs to albedo shower particle. Ignore for now.

    if (p == 0) return 0;
    
    // get original momentum from MCParticle
    
    float* pMC = p->GetMomentum();  
    float pTot = 0.;
    for ( int i=0; i < 3; i++) {
      pTot += pMC[i]*pMC[i];
    }
    pTot = sqrt(pTot);

    // convert to helical parameters

    float ptsqr = pMC[0]*pMC[0] + pMC[1]*pMC[1];
    float pt = sqrt(ptsqr);
    float ptInv = 1./pt;

    float cos_theta = pMC[2]/pTot;
    float abscth = fabs(cos_theta);
    float theta = acos(cos_theta);
    float lambda = TMath::Pi()/2. - theta;

    float phi = atan2(pMC[1],pMC[0]);
    if (phi < 0.) phi += 2.* TMath::Pi();

    McPart* parent = p->GetParnt();
    float* xMC = parent->GetPosition();
    float r = sqrt(xMC[0]*xMC[0]+xMC[1]*xMC[1]);

    int rc;

    // point at proper set of lookup tables

    LookUp2d* par[5];

    if ( abscth < m_parameters->getPolarInner() ) {
      for (int in=0; in < 5; in++) {
	par[in] = getBarrelTable(in);
      }
    }
    else if (abscth < m_parameters->getPolarOuter() ) {
      for (int o=0; o < 5; o++) {
	par[o] = getEndcapTable(o);
      }
    }

    // get resolution values from interpolation and smear parameters


    float ptInvDelta = par[ptI]->interpolateVal(abscth,pTot,&rc);
    float ptInvNew = ptInv + m_random->Gaus(0.,ptInvDelta);
    float ptNew = fabs(1./ptInvNew);

    float lambdaDelta = par[lambdaI]->interpolateVal(abscth,pTot,&rc);
    float lambdaNew = lambda + m_random->Gaus(0.,lambdaDelta);
    if (lambdaNew < -TMath::Pi()/2.) lambdaNew = -TMath::Pi()/2.;
    if (lambdaNew > TMath::Pi()/2.) lambdaNew = TMath::Pi()/2.;

    float rDelta = par[rI]->interpolateVal(abscth,pTot,&rc);
    float rNew = r + m_random->Gaus(0.,rDelta);

    float phiDelta = par[phiI]->interpolateVal(abscth,pTot,&rc);
    float phiNew = phi + m_random->Gaus(0.,phiDelta);
    if (phiNew < 0.) phiNew += 2.* TMath::Pi();

    float zDelta = par[zI]->interpolateVal(abscth,pTot,&rc);
    float zNew = xMC[2] + m_random->Gaus(0.,zDelta);

    // see if the charge sign has flipped

    float newCharge = (ptInv * ptInvNew < 0.) ? 
                        -p->GetCharge() : p->GetCharge();

    float thetaNew = TMath::Pi()/2. - lambdaNew;
    float pNew = ptNew/sin(thetaNew);


// get back to cartesian coordinates

    float newP[3], newX[3];
    
    newP[0] = pNew * sin(thetaNew) * cos(phiNew);
    newP[1] = pNew * sin(thetaNew) * sin(phiNew);
    newP[2] = pNew * cos(thetaNew);

    newX[0] = rNew * cos(phiNew);
    newX[1] = rNew * sin(phiNew);
    newX[2] = zNew;

    Track* tk = new Track(&newP[0],&newX[0],newCharge,index);

    return tk;

}

 SmearTrack::SmearTrack(GetParameters* gp){

  // constructor: sets up lookups tables, etc.

  m_random = new TRandom();
  m_parameters = gp;


  FILE* LTBFile;
  LTBFile = fopen(gp->getBarrelTableFile()->Data(),"r");


  GetTrackLookups* TkLookBarrel = new GetTrackLookups(LTBFile);
  TObjArray* barrelArray = TkLookBarrel->getLookupsArray();

  Int_t numPar = barrelArray->GetEntries();
  
  for (int i=0; i<numPar ; i++) {
    TString* parName = ((LookUp2d*)barrelArray->At(i))->getName();
    if ( parName->CompareTo("lambda",TString::kExact) == 0) {
      barrel_parNum[lambdaI] =((LookUp2d*)barrelArray->At(i));
    }
    else if ( parName->CompareTo("1/p_t",TString::kExact) == 0) {
      barrel_parNum[ptI]=((LookUp2d*)barrelArray->At(i));
    }
    else if ( parName->CompareTo("phi",TString::kExact) == 0) {
      barrel_parNum[phiI]=((LookUp2d*)barrelArray->At(i));
    }
    else if ( parName->CompareTo("z",TString::kExact) == 0) {
      barrel_parNum[zI]=((LookUp2d*)barrelArray->At(i));
    }
    else if ( parName->CompareTo("r",TString::kExact) == 0) {
      barrel_parNum[rI]=((LookUp2d*)barrelArray->At(i));
    }
    
  };

  printf("Barrel Smear parameters initialized n");
  
  if (!gp->getEndcapTableFile()->Contains("none")) {
    FILE* LTEFile;
    LTEFile = fopen(gp->getEndcapTableFile()->Data(),"r");
    
    GetTrackLookups* TkLookEndcap = new GetTrackLookups(LTEFile);
    TObjArray* endcapArray = TkLookEndcap->getLookupsArray();
    
    for (int j=0; j<numPar ; j++) {
      TString* parName = ((LookUp2d*)endcapArray->At(j))->getName();
      if ( parName->CompareTo("lambda",TString::kExact) == 0) {
	endcap_parNum[lambdaI]=((LookUp2d*)endcapArray->At(j));
      }
      else if ( parName->CompareTo("1/p_t",TString::kExact) == 0) {
	endcap_parNum[ptI]=((LookUp2d*)endcapArray->At(j));
      }
      else if ( parName->CompareTo("phi",TString::kExact) == 0) {
	endcap_parNum[phiI]=((LookUp2d*)endcapArray->At(j));
      }
      else if ( parName->CompareTo("z",TString::kExact) == 0) {
	endcap_parNum[zI]=((LookUp2d*)endcapArray->At(j));
      }
      else if ( parName->CompareTo("r",TString::kExact) == 0) {
	endcap_parNum[rI]=((LookUp2d*)endcapArray->At(j));
      }
      
    }
    printf("Endcap Smear parameters initialized n");
  }

};

int SmearTrack::lambdaI=0;
int SmearTrack::ptI = 1;
int SmearTrack::phiI =2;
int SmearTrack::zI = 3;
int SmearTrack::rI = 4;