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

//________________________________________
//  This is data header class
//  to store montecarlo simulation data
// 
// REMEMBER to call SetDefaults() after all parameter adjustment!
//          example is in script create_header.C
//
// To initialize fbias array call 
// ResetBias() or SetBias(0.0,0.0);
//
//
// Logs:
//
//  $Log: LDataHeader.cc,v $
//  Revision 1.24  2004/09/03 10:37:55  cvs
//  bias array is now possible to be loaded from parameters by external file
//
//  Revision 1.23  2004/08/25 08:05:00  cvs
//  there was a bug in deletion of fbias array. now corrected
//
//  Revision 1.22  2004/04/15 12:19:21  cvs
//  bugs...
//
//  Revision 1.21  2004/04/13 15:14:21  cvs
//  bugs...
//
//  Revision 1.20  2004/04/13 14:13:33  cvs
//  hopping is client, that can comunicate only throuth stdin/out
//
//  Revision 1.19  2004/04/09 13:21:05  cvs
//  added cputime logging...
//
//  Revision 1.18  2004/04/09 09:54:37  cvs
//  found why there are some peaks at bins e.g. 22,63,105,147... when TOF is drawn.
//  This was due to different number of bins in simulation and storage. To avoid this, both numbers of bins are set equal and less than MAX_BINS!
//
//  Revision 1.17  2004/02/09 14:26:25  cvs
//  added new LDataHeader version and force field also in streamer
//
//  Revision 1.16  2003/12/23 08:38:28  cvs
//  bug workaround in TLServer done.
//
//  Revision 1.15  2003/12/22 14:32:36  cvs
//  added simulation of laser -> LDataHeader version 5
//  modified hoppserver sigint interruption
//
//  Revision 1.14  2003/12/15 16:17:04  cvs
//  found a bug in SetDefaults. Time range was slightly miscalculated, so notches occured in all graphs...
//
//  Revision 1.13  2003/11/24 11:27:15  cvs
//  added version control
//
//  Revision 1.12  2003/11/11 17:33:03  cvs
//  LServer upgraded
//
//  Revision 1.11  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.10  2003/09/25 10:57:40  cvs
//  added force member to LDataHeader class
//
//  Revision 1.9  2003/08/28 14:41:18  cvs
//  corrected the bug in LData::Add(float,float,float,...). Description of bug: position index, that was obtained by simulation, was not shifted by -1 like was shifted the range of position in LDataHeader::SetDefaults.
//
//  Revision 1.8  2003/08/20 10:26:50  cvs
//  Added scripts directory with two utilities now:
//  modify_header which modify existing header and send it database
//  deliver_header program and script which deliver header to clients
//
//  libLComm was made ROOT compatible and GetBias(z) was added to LDataHeader
//
//  Revision 1.7  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.6  2003/07/15 14:19:35  cvs
//  corrected ranges when creating histograms. See LDataHeader::Defaults(). Added Draw(".bias.") to draw bias array
//
//  Revision 1.5  2003/07/09 13:27:39  cvs
//  Modified ranges for histograms stored in LData. Ranges are set by LDataHeader->SetDefaults() and MAX_BINS in Config.h
//
//  Revision 1.4  2003/06/26 14:46:16  cvs
//  found a bug in LSpace::GetEnergy() with fbias. Bug successfully resolved :)
//
//  Revision 1.3  2003/06/26 09:27:11  cvs
//  added log fields...
//

#include "LDataHeader.h"
#include <TROOT.h>//for gROOT
#include <TPad.h>//for gPad
#include <iostream>
#include <fstream>
#include <TRandom.h>

ClassImp(LDataHeader)

 LDataHeader::LDataHeader():TNamed("header","Monte Carlo simulation header"){
  force=-1;
  fbias=0;
  flasermode=0;
  fwspacex=fwspacey=fwspacez=7;//default values

  f_IsDefault=kFALSE;
}

 LDataHeader::LDataHeader(const char* filename):TNamed("header","Monte Carlo simulation header"){
  //constructor, which sets all data from a file. It uses Add(filename)
  fjobID=0;
  fbasex=fbasey=fbasez=1.0;
  falpha=0.0;
  fbias=0;
  flasermode=0;
  force=-1;//uninitalized
  fwspacex=fwspacey=fwspacez=7;//default values
  Add(filename);
  f_IsDefault=kFALSE;
}

 LDataHeader::LDataHeader(const int argc,const char** argv):TNamed("header","Monte Carlo simulation header"){
  //constructor, which sets all data from input line parameters, using the same names of parameters as in the file. uses Add(arvc,argv)
  fjobID=0;
  fbasex=fbasey=fbasez=1.0;
  falpha=0.0;
  fbias=0;
  flasermode=0;
  force=-1;//uninitialized
  fwspacex=fwspacey=fwspacez=7;//default values
  Add(argc,argv);
  f_IsDefault=kFALSE;
}

 LDataHeader::~LDataHeader(){
  //destructor
  //free memory used by fbias array
  if (fbias) {Tdelete(fbias,0);fbias=0;}
}

 void LDataHeader::Add(const char* filename){
  //sets all data from a file 
  FILE* fp=fopen(filename,"r");  
  if(!fp) {
    std::cerr<<"Configuracijska datoteka "<<filename<<" ne obstaja!"<<std::endl;
    return;
  }
  char str[256];
  while(fgets(str,256,fp)) {
    if(strncmp(str,"//",2) && strncmp(str,"#",1))
      Recognize(str);
  }
  fclose(fp);
};

 void LDataHeader::Add(const int argc,const char** argv){
  //sets all data from input line parameters, using the same names of parameters as in the file
  if(argc==1) {
    std::cerr<<"Error: No input line parametern";
    return;
  }
  for (int i=1;i<argc;i++) Recognize(argv[i]);  
};

 void LDataHeader::Recognize(const char* what){
  //Recognize parameters from input file
  //OBSOLETE use scripts create_header.C and modify_header.C instead!

  char label[256];
  char* value;
  
  strcpy(label,what);
  if (!index(label,'=')) return;
  value=&index(label,'=')[1];
  index(label,'=')[0]=0;
  //      printf("%s,,,%sn",label,value);
  
  TString labela(label);
  if (!labela.CompareTo("fjobid", TString::kIgnoreCase))
    {sscanf(value,"%ld",&fjobID);return;}
  if (!labela.CompareTo("fmax", TString::kIgnoreCase))
    {sscanf(value,"%d",&fmax);return;}
  if (!labela.CompareTo("fn1", TString::kIgnoreCase))
    {sscanf(value,"%d",&fN1);return;}
  if (!labela.CompareTo("force", TString::kIgnoreCase))
    {sscanf(value,"%d",&force);return;}
  if (!labela.CompareTo("fwspacex", TString::kIgnoreCase))
    {sscanf(value,"%d",&fwspacex);return;}
  if (!labela.CompareTo("fwspacey", TString::kIgnoreCase))
    {sscanf(value,"%d",&fwspacey);return;}
  if (!labela.CompareTo("fwspacez", TString::kIgnoreCase))
    {sscanf(value,"%d",&fwspacez);return;}
  if (!labela.CompareTo("fbasex", TString::kIgnoreCase))
    {sscanf(value,"%e",&fbasex);return;}
  if (!labela.CompareTo("fbasey", TString::kIgnoreCase))
    {sscanf(value,"%e",&fbasey);return;}
  if (!labela.CompareTo("fbasez", TString::kIgnoreCase))
    {sscanf(value,"%e",&fbasez);return;}
  if (!labela.CompareTo("fspacex", TString::kIgnoreCase))
    {sscanf(value,"%d",&fspacex);return;}
  if (!labela.CompareTo("fspacey", TString::kIgnoreCase))
    {sscanf(value,"%d",&fspacey);return;}
  if (!labela.CompareTo("fspacez", TString::kIgnoreCase))
    {sscanf(value,"%d",&fspacez);return;}
  if (!labela.CompareTo("fresolution", TString::kIgnoreCase))
    {sscanf(value,"%e",&fresolution);return;}
  if (!labela.CompareTo("ftemperature", TString::kIgnoreCase))
    {sscanf(value,"%e",&ftemperature);return;}
  if (!labela.CompareTo("fgamma", TString::kIgnoreCase))
    {sscanf(value,"%e",&fgamma);return;}
  if (!labela.CompareTo("fsigma", TString::kIgnoreCase))
    {sscanf(value,"%e",&fsigma);return;}
  if (!labela.CompareTo("falpha", TString::kIgnoreCase))
    {sscanf(value,"%e",&falpha);return;}
  if (!labela.CompareTo("flasermode", TString::kIgnoreCase))
    {sscanf(value,"%d",&flasermode);return;}
  if (!labela.CompareTo("flaserpar", TString::kIgnoreCase))
    {sscanf(value,"%e",&flaserpar);return;}
  Float_t temp;
  if (!labela.CompareTo("bias", TString::kIgnoreCase))
    {sscanf(value,"%e",&temp);SetBias(temp);return;}
  char filename[256];
  if (!labela.CompareTo("bias_arr", TString::kIgnoreCase))
    {sscanf(value,"%s",filename);//read bias array from file
    std::ifstream infile;
    infile.open(filename,ifstream::in
    if (!infile.is_open()){
      if (fverbose) std::cerr<<"Error: opening file:""<<filename<<""!n";
      return;}
    ResetBias();
    for (int i=0;i<=fspacez+1;i++){
      if (infile.good())
	infile >> fbias[i];
      else 
	if (fverbose) { std::cerr<<"Error: reading file:""<<filename<<""!n";
	return;}
    }
    infile.close();
    f_IsDefault=kFALSE;//fbias array has changed
    return;}
  
  if (fverbose) std::cerr<<"Info: unknown parameter:""<<what<<"" - Ignoring!n";
}

 void LDataHeader::PrintARG(){
  //print out string that is appropriate to initialize new header!
  std::cout<<"fjobid="<<fjobID;
  std::cout<<" fmax="<<fmax;
  std::cout<<" fn1="<<fN1;
  std::cout<<" force="<<force;
  std::cout<<" fwspacex="<<fwspacex;
  std::cout<<" fwspacey="<<fwspacey;
  std::cout<<" fwspacez="<<fwspacez;
  std::cout<<" fspacex="<<fspacex;
  std::cout<<" fspacey="<<fspacey;
  std::cout<<" fspacez="<<fspacez;
  std::cout<<" fbasex="<<fbasex;
  std::cout<<" fbasey="<<fbasey;
  std::cout<<" fbasez="<<fbasez;
  std::cout<<" fresolution="<<fresolution;
  std::cout<<" ftemperature="<<ftemperature;
  std::cout<<" fgamma="<<fgamma;
  std::cout<<" fsigma="<<fsigma;
  std::cout<<" falpha="<<falpha;
  std::cout<<" flasermode="<<flasermode;
  std::cout<<" flaserpar="<<flaserpar;
  //store bias array into text file
  const char* filename="bias_array";
  std::cout<<" bias_arr="<<filename<<std::endl;
  std::ofstream outfile;
  outfile.open(filename,ofstream::out
  if (!outfile.is_open()){
    if (fverbose) std::cerr<<"Error: opening file:""<<filename<<""!n";
    return;}
  for (int i=0;i<=fspacez+1;i++) outfile<<fbias[i]<<std::endl;
  outfile.close();
}

 void LDataHeader::CorrectRanges(Float_t& low,Float_t& up,Int_t n){
  //this is used to correct lower and upper range of the histogram
  //the problem is for eg. the coordinates of eight hopping sites 
  //are 1,2,3,4,5,6,7,8 If you want do draw histograms you have to
  //enter lower and upper range, which is 0.5 and 8.5
  //if you say 1 and 8 it is not correct. So this method correct the 
  //ranges, so that the center of histogram bin is the true coordinate
  Double_t delta=0.5*(up-low)/((Double_t)(n-1));
  low-=delta;
  up+=delta;
}
  
  
  
 void LDataHeader::SetDefaults(){
  //this set default parameters for data storing in LData
  //WARNING:call this method after having fbias filled!!!
  //time is measured from 0 
  //NEVER FORGET THIS!!!
  if (fmax>MAX_BINS){
    if (fverbose) std::cout<<"Info: fmax is too big! Reduced to "<<MAX_BINS<<"n";
    fmax=MAX_BINS;
  }  

  nbint=fmax;       
  rtmin=0;       
  rtmax=(fmax-1)*fresolution;
  CorrectRanges(rtmin,rtmax,nbint);

  //position in LDisorder is measured within 1 and fspacez
  //we shift it here, so the minimal value is 0
  if (fspacez>MAX_BINS) nbinp=MAX_BINS; else nbinp=fspacez;       
  rpmin=0;  
  rpmax=(fspacez-1)*fbasez;       
  CorrectRanges(rpmin,rpmax,nbinp);

  nbinv=nbinp;       
  Float_t maxjump=0.25*(fspacez*fbasez);//try
  Float_t maxvelocity=0.1*maxjump;
  rvmin=-maxvelocity;       
  rvmax=maxvelocity;       

  nbine=nbinp; 
  if (!fbias) {
    std::cerr<<"SetDefaults Error:: set bias or gradient before!n";
    return;
  }
  Float_t femax,femin;
  femax=femin=fbias[1];
  for (int i=2;i<=fspacez;i++) {
    if (femax>fbias[i]) femax=fbias[i];
    if (femin<fbias[i]) femin=fbias[i];
  }
  femax*=-1.0;
  femin*=-1.0;
  femax+=fsigma;
  femin-=fsigma;
  remin=femin;       
  remax=femax;       

  f_IsDefault=kTRUE;
}

 void LDataHeader::SetGradient(const Float_t gradient,const Float_t Uo){
  //calculate fbias array. It starts at the index, which is outside 
  //the hopping states, because the metal is at known bias, not organic!
  //U(z)=Uo-E*z    E...const
  //Uo=fbias[0] ... default is 0.0eV
  
  if (fspacez<1) {
    if (fverbose) std::cerr<<"Gradient error: use parameter SPACEZ before GRADIENT and BIAS!n";
    return;
  }
  if (!fbias) {
    Tnew(fbias,0,fspacez+1);
  }

  f_IsDefault=kFALSE;

  if (!fbasez>0.0) {
    if (fverbose) std::cerr<<"Gradient warning: use parameters SPACEZ,BASEZ before GRADIENT and BIAS!n";
    return;
  }
  
  Float_t dd=gradient*fbasez;
  Tclear(fbias,0,fspacez+1);
  fbias[0]=Uo;
  fbias[1]-=dd;
  for(Int_t i=2;i<=fspacez+1;i++)
    fbias[i]=fbias[i-1]-dd;
}

 TH1F* LDataHeader::GetBiasHist(){
  //see LDataHeader::SetDefults to se how the ranges are set
  Int_t nx = fspacez+2;
  Float_t min=-fbasez;
  Float_t max=fbasez*fspacez;
  
  CorrectRanges(min,max,nx);
  
  const char* pname="bias";    
  TH1F *h1=0;
  //check if a hist with identical name exist
  TObject *h1obj = gROOT->FindObject(pname);
  if (h1obj && h1obj->InheritsFrom("TH1F")) {
    h1 = (TH1F*)h1obj;
    h1->Reset();
  }
  
  if (!h1) {
    h1 = new TH1F(pname,"Bias vs Position",nx,min,max);
    h1->SetXTitle("Position [Ang.]");
    h1->SetYTitle("Bias Potential [eV]");  
  }

  // Fill the histogram
  for (Int_t binx =1;binx<=nx;binx++) h1->SetBinContent(binx,fbias[binx-1]);
  return h1;
}


 void LDataHeader::SetBias(const Float_t bias,const Float_t Uo){
  //calculate and correct fbias array. It starts at the z=0.
  //it is advised to adjust fbias[0] to contain the correct 
  //value
  //U(z)=Uo-E*z    E...const
  //Uo=fbias[0] ... default is 0.0eV
  
  Float_t d=fbasez*(Float_t)(fspacez+1);
  if (d>0.0) SetGradient(-bias/d,Uo);
  else SetGradient(0.0,Uo);
}
      

 void LDataHeader::SetBias(Float_t* bias){
  //copy array to new one
  ResetBias();//used to reset
  for(int i=0;i<=fspacez+1;i++)
    fbias[i]=bias[i];
}

 Float_t  LDataHeader::GetBias(const Int_t z){
  //return bias potential for sites
  //argument z is site index [0..fspacez+1]
  if (z<1) return fbias[0];
  if (z>fspacez) return fbias[fspacez+1];
  return fbias[z];
}

 Float_t  LDataHeader::GetBias(const Float_t z){
  //return bias potential in between points
  //argument z is space coordinate not site index [-fbasez..fspacez*fbasez]
  //see Defaults() for such ranges

  Float_t zz=z+fbasez;//just to simplify calculations
  if (zz<=0.0) return fbias[0];
  if (zz>=((fspacez+1)*fbasez)) return fbias[fspacez+1];
  Float_t dd=(zz/fbasez);
  Int_t i=(Int_t) (zz/fbasez);
  dd-=(Float_t) i;
  return fbias[i]+(fbias[i+1]-fbias[i])*dd;
}

void 
 LDataHeader::SetLaserMode(const Int_t mode,const Float_t par){
  //set the mode to simulate laser operation. Available modes are:
  // 0 - no laser simulation, all are generated at time zero
  // 1 - gaus laser simulation, parameter [0] is sigma of gaus pulse duration
  //     sigma is in units of fmax, floating points are allowable
  flasermode=mode;
  if (mode==1){
    if (!par) std::cerr<<"Error in SetLaserMode(): No sigma set! Ignoring!n";
    else{ 
      if (par>0.0) flaserpar=par*fresolution;
      else std::cerr<<"Error in SetLaserMode(): Sigma negative! Ignoring!n";
    }
  }
}

Float_t 
 LDataHeader::CalcIndexToTime(const Int_t idx) const{
  //this calculate real time from LDisorder experiment
  //here is the laser simulation
  Float_t offset=0.0;
  if (flasermode==1) offset=gRandom->Gaus(flaserpar*3,flaserpar);
  return (Float_t (idx-1))*fresolution+offset;
}

void 
 LDataHeader::Browse(TBrowser*){
  Inspect();
  gPad->Update();
}

void 
 LDataHeader::Streamer(TBuffer &R__b){
  //special streamer for LDataHeader
  //pazi fbias array
  if(R__b.IsReading()){
    Version_t R__v;
    R__v=R__b.ReadVersion();
    TNamed::Streamer(R__b);
    R__b>>fjobID;
    //    std::cout<<"received jobID:"<<fjobID<<std::endl;
    R__b>>fmax;
    R__b>>fN1;
    if (R__v > 5) {
      R__b>>force;
    }
    if (R__v > 3) {
      R__b>>fwspacex;
      R__b>>fwspacey;
      R__b>>fwspacez;
    }
    R__b>>fspacex;
    R__b>>fspacey;
    R__b>>fspacez;
    R__b>>fbasex;
    R__b>>fbasey;
    R__b>>fbasez;
    R__b>>fresolution;
    R__b>>ftemperature;
    R__b>>fgamma;
    R__b>>fsigma;
    R__b>>falpha;
    if (R__v > 4) {
      R__b>>flasermode;
      R__b>>flaserpar;
    }
    R__b>>rtmin;
    R__b>>rtmax;
    R__b>>rpmin;
    R__b>>rpmax;
    R__b>>rvmin;
    R__b>>rvmax;
    R__b>>remin;
    R__b>>remax;
    R__b>>nbint;
    R__b>>nbinp;
    R__b>>nbinv;
    R__b>>nbine;  

    if (fspacez){
      if (fbias) Tdelete(fbias,0);
      else Tnew(fbias,0,fspacez+1);
      R__b.ReadFastArray(fbias,fspacez+2);
    }
    
    if (R__v > 6) {
      R__b>>f_IsDefault;
    } else f_IsDefault=kTRUE;

  } else {
    Version_t R__v=LDataHeader::IsA()->GetClassVersion();
    R__b.WriteVersion(LDataHeader::IsA());
    TNamed::Streamer(R__b);
    //    std::cout<<"send jobID:"<<fjobID<<std::endl;
    R__b<<fjobID;
    R__b<<fmax;
    R__b<<fN1;
    if (R__v > 5) {
      R__b<<force;
    }
    if (R__v > 3) {
      R__b<<fwspacex;
      R__b<<fwspacey;
      R__b<<fwspacez;
    }
    R__b<<fspacex;
    R__b<<fspacey;
    R__b<<fspacez;
    R__b<<fbasex;
    R__b<<fbasey;
    R__b<<fbasez;
    R__b<<fresolution;
    R__b<<ftemperature;
    R__b<<fgamma;
    R__b<<fsigma;
    R__b<<falpha;
    if (R__v > 3) {
      R__b<<flasermode;
      R__b<<flaserpar;
    } 
    R__b<<rtmin;
    R__b<<rtmax;
    R__b<<rpmin;
    R__b<<rpmax;
    R__b<<rvmin;
    R__b<<rvmax;
    R__b<<remin;
    R__b<<remax;
    R__b<<nbint;
    R__b<<nbinp;
    R__b<<nbinv;
    R__b<<nbine;  

    if (!fbias){
      Tnew(fbias,0,fspacez+1);
      ResetBias();
    }
    R__b.WriteFastArray(fbias,fspacez+2);

    if (R__v > 6) {
      if (!f_IsDefault)
	std::cerr<<"LDataHeader::Streamer: Warning: Have you called SetDefaults()?  or have you forgot to set Bias?"<<std::endl;
      R__b<<f_IsDefault;
    } 
  }
}