// Author: Egon Pavlica <http://www.p-ng.si/~pavlica/>

//________________________________________
//  This is IV class
//  to simulate IV caracteristics
//  using monte carlo simulation
//  stored in LData
// 
//  the scenario is the follows:
//  1. get LData and create this class with
//  2. set input parameters with Set*
//  3. run with Adjust bias
//  4. get new header and do MC simulations
//  5. delete this class to release LData
//
//  See the documentation in doc directory for more info...
//
// Logs:
//
//  $Log: LIV.cc,v $
//  Revision 1.7  2003/11/03 13:54:40  cvs
//  TLServer is now threaded, LData and LDataHeader have now right Draw and Paint methods, LIV is controlled and comments were added.
//
//  Revision 1.6  2003/08/30 11:28:52  cvs
//  simplified and checked LIV calculations
//
//  Revision 1.5  2003/07/25 11:36:26  cvs
//  simulation of IV measurements is done with iv program
//  some minority fixes to LDisorder - corrected ft to match the recorded position
//
//  Revision 1.4  2003/07/15 15:20:00  cvs
//  started LIV development, including doc file.
//
//  Revision 1.2  2003/06/26 09:27:11  cvs
//  added log fields...
//

#include "LIV.h"
#include <iostream>
#include <TProfile.h>
#include <TH1F.h>
ClassImp(LIV)

 LIV::LIV(LData* data){
  //input simulation data is used to calculate charge distribution 
  //which should be present in constant current.
  if (!data)
    std::cerr<<"Error:LIV no LData classn"<<std::endl;  
  fdata=data;
  //create new header with different fbias array
  fheader=(LDataHeader*)fdata->GetHeader()->Clone();
  fheader->SetJobID(0);//this will be new job  
  //fheader->Dump();

  fN=fheader->fspacez+2;//n of bins in fheader->fbias 

  //calculation arrays
  Tnew(ffield,0,fN);  Tclear(ffield,0,fN);
  Tnew(fbias,0,fN);  Tclear(fbias,0,fN);
  Tnew(fcharge,1,fN-1); Tclear(fcharge,1,fN-1);
  Tnew(fv,1,fN-1); Tclear(fv,1,fN-1);
  Tnew(fz,1,fN-1); Tclear(fz,1,fN-1);

  fe0=kEPSILON0;//maybe is in e/VAng. 

  //posible memory leak here!!!
  TProfile* his=(TProfile*) data->GetHistogram(".avp.");
  TH1F* bias=(TH1F*) data->GetHistogram(".bias.");

  fdz=bias->GetBinWidth(1);//hope they are all the same width

  for (int i = 1;i<fN;i++){
    fz[i]=bias->GetBinCenter(i);
    fv[i]=his->GetBinContent(his->GetXaxis()->FindBin(fz[i]));
  }

}

 LIV::~LIV(){
  Tdelete(fz,1);
  Tdelete(fv,1);
  Tdelete(ffield,0);
  Tdelete(fbias,0);
  Tdelete(fcharge,1);
}


Float_t 
 LIV::Chi2(const Float_t* a,const Float_t* b,const Int_t n)const {
  //calculate chi2 between two arrays 1-base-index
  //return sqrt(chi2)
  Double_t sum=0.0;
  Double_t del;
  for(int i=1;i<=n;i++){
    del=(Double_t)(a[i]-b[i]);
    sum+=del*del;
  }
  return (sqrt(sum));
}


Float_t
 LIV::AdjustBias(){
  //steps:
  // 1-charge density is calculated to generate
  //   constant current - here is possible to make some intelligent
  //   method to calculate good charge density. This used is very simple
  // 2-electric field is calculated from charge density
  // 3-fbias is calculated and adjusted to match bias
  // v(x)=vo.v'(x)   vo is returned
  // v'(x) is obtained from LData

  if (fI==0.0) std::cerr<<"WARNING! LIV::AdjustBias(): use SetI() to set currentn";
  if (fe==0.0) std::cerr<<"WARNING! LIV::AdjustBias(): use Sete() to set en";
  if (fS==0.0) std::cerr<<"WARNING! LIV::AdjustBias(): use SetS() to set Sn";

  fI/=1.6e-19; //if curent is entered in A

  for(Int_t i=1;i<fN;i++)
    if (fv[i]>0.0) fcharge[i]=1.0/fv[i];


  //  ffiled[0] is set elsewhere
  //  ffiled[fN] is set elsewhere
  ffield[1]=fcharge[1];
  for(Int_t i=2;i<fN;i++)
    ffield[i]=ffield[i-1]+fcharge[i];

  // fbias[0] is set elsewhere
  // fbias[fN] is set elsewhere
  fbias[1]=ffield[1];
  for(Int_t i=2;i<fN;i++)
    fbias[i]=fbias[i-1]+ffield[i];

  //we can now calc v0;
  Float_t dd;//those long correction to fbias
  Float_t kk;
  std::cout<<"1-fv0 should be 0:"<<fvo<<"n";
  dd=ffield[0]+ffield[fN];
  std::cout<<"2-this should be Eo+En:"<<dd<<"n";
  std::cout<<"3-dz should not be 0:"<<fdz<<"n";
  dd*=fdz;
  dd=fbias[fN]-fbias[0]+dd;
  std::cout<<"4-this should be Un-Uo+(Eo+En)dz:"<<dd<<"n";
  if (dd==0.0) std::cerr<<"WARNING! LIV::AdjustBias(): use SetUn() or others to set potential differencen";
  dd/=fbias[fN-1];//note here that fN=N+1 regarding documentation
  std::cout<<"5-this should be prev. value divided by "<<fbias[fN-1]<<" :"<<dd<<"n";
  kk=fI*fdz*fdz;
  std::cout<<"6-I*dz**2/eo:"<<kk<<"n";
  kk/=(fS*fe*fe0);
  std::cout<<"7-previous divided by S*e*eo="<<(fS*fe*fe0)<<" is :"<<kk<<"n";
  fvo=-1.0*kk/dd;
  std::cout<<"8-fvo is -[7]/[5]:"<<fvo<<"n";
  //here we adjust fbias by multiplying
  for(Int_t i=1;i<fN;i++) fbias[i]*=dd;
  //here we adjust fbias by shifting
  Float_t shift=fbias[0]-ffield[0]*fdz;
  for(Int_t i=1;i<fN;i++) fbias[i]+=shift;
  
  //prepare new header
  fheader->SetBias(fbias);

  //calculate chi2
  return Chi2(fheader->GetBias(),fdata->GetHeader()->GetBias(),fN);
  
}


 TH1F* LIV::GetTOF(){
  //returns the TOF plot in real time units, obtained from LIV model
  if (fvo==0.0) return 0;
  float to=Getto();
  TH1F* his=(TH1F*)fdata->GetHistogram(".tof.");
  int n=his->GetXaxis()->GetNbins();
  float start=his->GetXaxis()->GetBinLowEdge(1);
  float end=his->GetXaxis()->GetBinLowEdge(n)+his->GetXaxis()->GetBinWidth(n);
  his->GetXaxis()->Set(n,start*to,end*to);
  his->GetXaxis()->SetTitle("Time [s]");
  return his;
}