#include "SiRDPScen.h"
#define NMAX 10000
ClassImp(SiRDPScen)

// Double_t SimpleFlux(Double_t *x,Double_t *par)
// {
//   Int_t NumPar=(Int_t) par[0];
  

//   return (x[0]+par[0]);
// }

 Int_t SiRDPScen::ReadTrapsFile(Char_t *trapfile_name)
{ 
  //Function reads traps file.
  //Each trap in the trap file contains one line per trap e.g.:
  //A 1 0.02 0.000 0.1 1.11 293 -1. 0. 0. 0
  //The parameters describing the trap in the file are:
  //1. the trap name in this case : A
  //2. the dynamics (1st , 2nd)   : 1
  //3. generation rate [1/cm]     : 0.02
  //4. generation rate [1/cm]     : 0.000 
  //           in case of second order process reaction A+B->C
  //           is parameter 3. generation rate of A and 4. of B
  //5. reaction constant k        : 0.1
  //           in case of 1st order process: k=1/tau
  //6. activation energy          : 1.11
  //7. reference temperature      : 293
  //8. contribution to the space charge of A: -1.
  //9. contribution to the space charge of B: 0.
  //10. contribution to the space charge of C: 0.
  //11. is defect sensitive to bias: 0
  FILE *in;
  Int_t ok;
  SiRDPTrap temp;
  if( (in=fopen(trapfile_name,"r+"))==NULL)
    {
      printf("Can not open the file: %s",trapfile_name);
      return 1;
    }
  while(!feof(in))
    {
      ok=temp.ReadTrapFile(in);
        if(!ok) AddTrap(&temp);
    }
}

 void SiRDPScen::Link(Int_t s1,Int_t s2,Int_t d1,Int_t d2)
{
  //The function links different defects, so that one can describe 
      Int_t NumS;
      TArrayF *Ns,*Nd;
      
      switch(s2)
	{
	case 0:
	  Ns=&(GetTrap(s1)->Na); 
	  NumS=(GetTrap(s1)->Na).GetSize(); break;
	case 1:
  	  Ns=&(GetTrap(s1)->Nb); 
	  NumS=(GetTrap(s1)->Nb).GetSize(); break;
	case 2:
  	  Ns=&(GetTrap(s1)->Nc);
	  NumS=(GetTrap(s1)->Nc).GetSize(); break;
	}

      switch(d2)
	{
	case 0:
	  Nd=&(GetTrap(d1)->Na); break;
	case 1:
  	  Nd=&(GetTrap(d1)->Nb); break;
	case 2:
  	  Nd=&(GetTrap(d1)->Nc); break;
	}

      Nd->Adopt(NumS,Ns->GetArray());
}


 SiRDPScen::SiRDPScen(Char_t *scenario_name)
{
  Nstep=0;
  Niter=0;
  if(scenario_name!=0)
  ReadScenario(scenario_name,1);
  Traps=new TClonesArray("SiRDPTrap",10);
  Traps->BypassStreamer(kFALSE); 
  Fluence.Set(Niter+1);
  Time.Set(Niter+1);
}

 Int_t SiRDPScen::AddTrapLinks(Int_t S, Int_t Sd,Int_t D,Int_t Wd)
{
Int_t num=List.GetSize();
List.Set(num+4);
 if(S>GetEntries() || S<0) {printf("Index error!n"); return 1;}
 if(D>GetEntries() || D<0) {printf("Index error!n"); return 2;}
 if(Sd>2 || Sd<0 || Wd>2 || Wd<0) 
   {printf("Index error!n"); return 3;}
 List[num]=S;
 List[num+1]=Sd;
 List[num+2]=D;
 List[num+3]=Wd;
 Link(S,Sd,D,Wd); 
 //Link(0,1,1,1); 
 return 0;
}

 void SiRDPScen::Delete()
{
 Int_t num=GetEntries();
 Int_t i=0;
 
    Traps->Delete();      

  Fluence.Reset();
}

 void SiRDPScen::Reset()
{
 Int_t num=GetEntries();
 Int_t i=0;
 
  for(i=0;i<num;i++)
    ((SiRDPTrap *)Traps->At(i))->Reset();      

  Fluence.Reset();
}

 void SiRDPScen::ShowLinks()
{
  Int_t i=0,j=0,k=0;
  TString t1,t2;
  Char_t n1,n2;
  Int_t num=List.GetSize();
  for(i=0;i<num;i+=4)
    {
      t1=GetTrap(List[i])->GetName();
      t2=GetTrap(List[i+2])->GetName();
      
      switch(List[i+1])
	{
	case 0: n1='a'; break;
        case 1: n1='b'; break;
	case 2: n1='c'; break;
	}
      switch(List[i+3])
	{
	case 0: n2='a'; break;
	case 1: n2='b'; break;
	case 2: n2='c'; break;
	}


       printf("Trap Link: %s(%d)%c -> %s(%d)%cn",(const char *) t1,List[i],n1,(const char *) t2,List[i+2],n2);
	
	k++;
    }
}


 void SiRDPScen::AddTrap(SiRDPTrap *trap)
{ 
  SiRDPTrap *tmp;
  Int_t last=(Int_t)Traps->GetLast();
  Traps->New( last+1 );
  tmp=(SiRDPTrap *)Traps->At( last+1 ); 
  trap->Copy(tmp);
  if(Nstep!=0) 
    {
      //printf("Warning: the content of Nx and Ny field will be deletedn");
           tmp->Na.Set(Niter+2); tmp->Na.Reset(-1); tmp->Na[0]=0;
           tmp->Nb.Set(Niter+2); tmp->Nb.Reset(-1); tmp->Nb[0]=0;
           if(tmp->GetModel()==2) 
	     {tmp->Nc.Set(Niter+2); tmp->Nc.Reset(-1); tmp->Nc[0]=0;}	
    }
}

 void SiRDPScen::ListTraps(Bool_t short_desc)
{ 
  Int_t i;
  for(i=0;i<=Traps->GetLast();i++)
    {
      printf("TRAP SERIAL NUMBER : %dn",i);
      if(short_desc)
      printf("Trap name : %sn",(const char *)((SiRDPTrap *) Traps->At(i))->GetName());
	else      
      ((SiRDPTrap *) Traps->At(i))->Info();
      printf("*******************************n");
    }

}


 SiRDPScen::~SiRDPScen()
{
  if(Traps!=NULL) delete Traps;
}

 Int_t SiRDPScen::ReadScenario(Char_t *file,Bool_t print)
{
  FILE *in;
  Int_t ok=1;
  //Opening the scenario file for reading 
  if((in=fopen(file,"r+"))==NULL) 
    {printf("Cannot open the scenario file!n"); return 1;}
  
//allocating the default storage

  et.Set(NMAX);
  flu.Set(NMAX);
  temp.Set(NMAX); 
  bias.Set(NMAX); 
  step.Set(NMAX);

//reading the data and checking its integrity

  while ((ok=fscanf(in,"%e %e %f %f %fn",&et[Nstep],&flu[Nstep],&temp[Nstep],&bias[Nstep],&step[Nstep]))!=EOF) 
   {
if(print) 
     printf("time int.=%8.3e : Fluence=%8.3e [1/cm2] : T=%8.3f [K] : U=%8.3f [V] : step=%8.3f n", et[Nstep],flu[Nstep],temp[Nstep],bias[Nstep],step[Nstep]);

     if(et[Nstep]<0) 
       { 
	 printf("Step %d: Check time interval! Error! n",Nstep); 
         Nstep=0;
	 return 2;
       }
     if(flu[Nstep]<0) 
       {
	 printf("Step %d: Check fluence! Error! n",Nstep); 
	 Nstep=0;
	 return 3;
       }

     if(temp[Nstep]<0 || temp[Nstep]>1500) 
       {
	 printf("Step %d: Check temperature! Error! n",Nstep); 
	 Nstep=0;
	 return 4;
       }
     if(bias[Nstep]<0 || bias[Nstep]>1000) 
       {
	 printf("Step %d: Check bias! Error! n",Nstep); 
	 Nstep=0;
	 return 5;
       }
     if(step[Nstep]<0 || step[Nstep]>et[Nstep]) 
       {
	 printf("Step %d: Check time step! Error! n",Nstep); 
	 Nstep=0;
	 return 6;
       }

   
	   
     Niter+=(Int_t)(et[Nstep]/step[Nstep]);

     Nstep++;
   }

  et.Set(Nstep);
  flu.Set(Nstep);
  temp.Set(Nstep); 
  bias.Set(Nstep); 
  step.Set(Nstep);
fclose(in);
return(0);
}   

 void  SiRDPScen::DoStep()
{
  Int_t num=GetEntries();
  SiRDPTrap *tmp;
  Int_t i,j;  
 
  for(i=0;i<num;i++)
    ((SiRDPTrap *)Traps->At(i))->DoStep();      

}



 void  SiRDPScen::Set(Int_t spos)
{
  Int_t num=GetEntries(); 
  SiRDPTrap *tmp;
  for(Int_t i=0;i<num;i++)
    {
      tmp=((SiRDPTrap *)Traps->At(i));
      tmp->SetFlux(flu[spos]/et[spos]);
      tmp->SetU(bias[spos]);
      tmp->SetT(temp[spos]);
      tmp->SetTimeStep(step[spos]);
    }
}


 void SiRDPScen::DoScenario(Int_t print)
{
  Int_t i=0,k=0;
  Int_t spos=0;
  Int_t siter;
  Float_t flux=0,timestep;
  Set(spos);
  //  ListTraps();
  siter=(Int_t) (et[spos]/step[spos]);
  flux=flu[spos]/et[spos];
  timestep=step[spos];
  if(print) printf("Looping over the scenario::::: n");
  if(print) printf("STEP 0 (Iter=%4.0d):",siter);
  for(i=0;i<Niter;i++)
    {
      if(k>siter) 
	 { 
	   spos++; 
	   Set(spos);
           flux=flu[spos]/et[spos];
	   timestep=step[spos];
	   siter=et[spos]/step[spos];
	   if(print) printf("nSTEP %d (Iter=%4.0d):",spos,siter);	   
	   k=0;
	 } 
   if(print)   if(i%((Int_t)(siter/20))==0) printf(".");
      DoStep();
      if(i>0) 
	{
	  Fluence[i]=Fluence[i-1]+flux*timestep;
	  Time[i]=Time[i-1]+timestep;
	}
      k++;
    }
   if(print) printf("nEnd of scenario!n");

}