//	$Log: BrDigitizeTof.cxx,v $
//	Revision 1.2  1998/03/06 22:10:01  videbaek
//	Working update
//
//	Revision 1.1.1.1  1998/03/04 21:33:37  brahmlib
//	Brat tof
//
//	$Id: BrDigitizeTof.cxx,v 1.2 1998/03/06 22:10:01 videbaek Exp $
//
//   File BrDigitizeTof.cxx
//
////////////////////////////////////////////////////////////////
//
//	BrDigitizeTof
//
//      Digitisation Module class for time-of flight detectors
//      which have each slat instrumented with two tubes at each end.
//      The module expect to get geant hits as input and produce digitized
//      Tof data BrDigTof. These could also have been labeled RawTof. We
//      may want to distuinguish between digitiaed as used in Brat and
//      raw dat as written by DAQ ?
// 
//
////////////////////////////////////////////////////////////////
//
#include "math.h"
#include <stdlib.h>
#ifdef UNIX
#include <unistd.h>
#endif
#include <iostream.h>

//
//  Root includes
//


#ifndef ROOT_TMath
  #include "TMath.h"
#endif
#include "TSonataMath.h"
#ifndef ROOT_TRandom
  #include "TRandom.h"
#endif

#include "BrDigitizeTof.h"
#include "BrDetectorParamsTof.h"
#include "TSonataMath.h"
#include "BrDetectorVolume.h"
#include "BrEventNode.h"
#include "BrDataTable.h"
#include "BrGeantHit.h"
#include "BrDigTof.h"

ClassImp(BrDigitizeTof)


 BrDigitizeTof::BrDigitizeTof()
{
// Constructor. 
}

 BrDigitizeTof::BrDigitizeTof(Text_t *Name,Char_t *Title)
  :BrModule(Name,Title)
{
fParamsTof_p = NULL;
}

 void BrDigitizeTof::SetDetectorParamsTof(BrDetectorParamsTof* par)
{
  // Set the Tof parameters with the values of the external parameters
  // referenced by pointer
  // Note that no owner ship is transferred. A copy is made.
  // 
  // 
  if( fParamsTof_p == NULL)
    fParamsTof_p = new BrDetectorParamsTof("Dummy","TofINT");
  *fParamsTof_p = *par;
  
}

 void BrDigitizeTof::SetDetectorParamsTof(const BrDetectorParamsTof& par)
{
  // Set the Tof parameters with the values of the external parameters.
  // Note that no owner ship is transferred. A copy is made.
  // There is no real need for the copy code here. Both Name and Title
  // will be copied (it seems)     
  if( fParamsTof_p == NULL)
    fParamsTof_p = new BrDetectorParamsTof("Dummy","ParamsTof");
  *fParamsTof_p = par;
  
}



 BrDigitizeTof::~BrDigitizeTof(){
  if( fParamsTof_p != NULL)
    delete fParamsTof_p;
};


 void BrDigitizeTof::Init(){
  // Initialize entry point. This is not called at creation time but
  // This would be a good place to setup histograms, statistics variables etc.
  // This could also be done at the time of the constructor. This might likely be more
  // Usefull since one has the opportunity to setup varibales before this call
  // (e.g. detector parameters or alike).
}

 void BrDigitizeTof::Event(BrEventNode* InputTable, BrEventNode* OutputTable)
{

  // Event method to be called once per event. 
  // This performes the actual slow simulation
  // of the hits on the Tof.
  //
  // Description of method:
  // Each slat can in fact be hit by multiple hits
  // 1) First a loop over all hits is made to see how many
  //    hits there are per slat. If the is only a single (charged) hit
  //	  The adc and tdc value is calculated from the position, de/dx and flight time
  // 2) For multiple hits each end of the tube is dealt with seperately.
  //	  The fastest signals is found (at each end). If the following hits are within
  //    the resolving time / gate width (100n sec) the adc values are
  //    added to thos from the first hit.
  //
 
  BrDetectorParamsTof *par;
  BrDetectorVolume    *vol;
  
  Int_t NumHits;
  BrDataTable* GeantHits;
  BrGeantHit* ghit_p;
  BrDigTof *digtof_p;
  if(DebugLevel()>0){ 
    printf("Inside Digitize of %sn",GetName());;
    InputTable->ListObjects();
  }
  //
  // Get volume information
  //
  vol = GetDetectorVolume();
  if(vol == NULL){
    cout << "Volume must be set " << endl;
    return ;
  }
  //
  // Get simulation parameters
  //
  if(DebugLevel()>0 ) vol->ListParameters();
  par = GetDetectorParamsTof();
  if(par == NULL){
    cout << "Parameters must be set " << endl;
    return ;
  }
  if( DebugLevel() > 2){
    cout << " Getting Parameters" << endl;
    par->ListParameters();
    cout << "here" << endl;
  }
  //
  // Decode the inputtable
  // The tof hits come from a table with e.g. the name "GeantHits H1"
  //
  Char_t TableName[32];
  sprintf(TableName,"GeantHits %sx0", GetName());
  GeantHits = (BrDataTable*)InputTable->GetObject(TableName);
  
  if(GeantHits == 0 ) {
    if(DebugLevel() > 0){
      cout << "No hits for " << GetName() << endl;
    }
    return;
}
  
  NumHits = GeantHits->Entries();
  //
  // ObjectList returns TObjectArray
  // Prepare the output table
  //
  sprintf(TableName,"BrDigTof %sx0",GetName());
  //
  BrDataTable *DigitizedTof = new BrDataTable(TableName);
  
  
  OutputTable->AddObject(DigitizedTof);
  
  int noslats = fParamsTof_p->GetNoSlats();
  
  int ip;     //
  int slat;
  int nohits; // No of hits for a given slat
  int ih;     // Local hit identifier.
  
  int* slattable = new int [noslats];
  BrGeantHit** hittable  = new BrGeantHit* [NumHits];
  for(ip=0;ip<noslats;ip++) slattable[ip]=0;


  if(DebugLevel()>1){
    cout << "No of Input hits = " << NumHits << endl;
  }
  //
  // Loop over all Slats;
  // Check on hits and create TofDig according to the algorithm given
  // in the calss description.
  //
  BrGeantHit* hit_p;
  int hitno = 0, nlookup, n;
  Float_t adc, tdc, t0, ypos, length, dl, dedx, decaylength;
  Float_t a,b, time, vspeed, tup, tdown, adc_up, adc_down;
  for(int is=0; is< noslats; is++){
    if(DebugLevel()>1){
      cout << "Check Slat No. " << is << endl;
    }
    //
    // Find all hits for this slat
    // and fill the pointer table.
    //

    for(ih=0; ih < NumHits; ih++){
      ghit_p = (BrGeantHit*) GeantHits->At(ih);
      //
      // this kind of check should really be done using assert(ghit_p)
      //
      if(ghit_p == NULL){
	cout <<  "**error ought not to happen" << endl;
	break;
      }
      slat = ghit_p->Isub()-1;
      if(DebugLevel()>2){
	cout <<  "Slatno " << slat << " de/dx " << ghit_p->Dedx() << endl;
      }
      if (slat == is){
	if(ghit_p->Dedx() > 0.0){
	  if(DebugLevel()>1){
	    cout <<  "Add Slatno " << slat << endl;
	  }
	  slattable[slat]++;
	  hit_p = ghit_p;
	  hittable[hitno] = hit_p;
	  hitno++; 
	} 
      }
    }      
    //
    // Single hit
    //
    if(DebugLevel()>3){
      cout << "Slattable " << slattable[is] << endl;
    }
    if(slattable[is] == 1){
      digtof_p = new BrDigTof();
      DigitizedTof->Add(digtof_p);
      digtof_p->SetSlatno(is+1);

      hit_p = hittable[0];

      length = fParamsTof_p->GetScintLength();
      dedx = hit_p->Dedx();
      ypos = hit_p->LocalYPosIn();
      decaylength = fParamsTof_p->GetAttLength();
      dl   = TSonataMath::Exp(-(length-ypos)/decaylength);
      adc  = dedx * fParamsTof_p->GetADCGainUp()
	     / fParamsTof_p->GetMeVperMIP()*0.001*dl;
      //
      // It would be nice to to refer to a global variable within root
      // but to a wrapped function in a specific namespace.
      //
      gRandom->Rannor(a,b);
      adc = adc*(1.0 + a * fParamsTof_p->GetSigmaAdc());
      digtof_p->SetAdcUp(adc);
      dl   = TSonataMath::Exp(-(length+ypos)/decaylength);
      adc  = dedx * fParamsTof_p->GetADCGainDown()
	     / fParamsTof_p->GetMeVperMIP()*0.001*dl;
      adc = adc*(1.0 + b * fParamsTof_p->GetSigmaAdc());
      digtof_p->SetAdcDown(adc);

      gRandom->Rannor(a,b);
      t0  = hit_p->Tof() + (length - ypos) * fParamsTof_p->GetSignalSpeed();
      t0 += a* fParamsTof_p->GetSigmaTime();
      time= t0/fParamsTof_p->GetTdcConv()
	+fParamsTof_p->GetTdcOffset();
      digtof_p->SetTdcUp(time);

      t0  = hit_p->Tof() + (length + ypos) * fParamsTof_p->GetSignalSpeed();
      t0 += a* fParamsTof_p->GetSigmaTime();
      time= t0/fParamsTof_p->GetTdcConv()
	+fParamsTof_p->GetTdcOffset();

      digtof_p->SetTdcDown(time);
    }
    //
    //  Multiple hits on same slat
    //
    else if(slattable[is] > 1){
      length = fParamsTof_p->GetScintLength();
      decaylength = fParamsTof_p->GetAttLength();
      nlookup = slattable[is];
      vspeed  = fParamsTof_p->GetSignalSpeed();
      //
      // Find the first time for the up and down Tube.
      //
      tup = tdown = 9999.0;
      for(n=0;n<nlookup;n++){
	    hit_p = hittable[n];
	    ypos = hit_p->LocalYPosIn();
	    time  = hit_p->Tof() +(length - ypos) * vspeed;
		cout << n << " " << hit_p->Tof() << " " << time << " "
			<< hit_p->Dedx() << " " << ypos << endl;
	    if(time < tup)
	      tup = time;
	    time  = hit_p->Tof() +(length + ypos) * vspeed;
	    if(time < tdown)
	      tdown = time;
      }
      //
      // Add the new BrDigTof to data.
      // For the time the fastest signal at each end is used.
	  // This code does not simulate any slewing behaviour
	  // For  adc signals the total is added up for time's within 100nsec
	  // of the fastest time.
	  //
      digtof_p = new BrDigTof();
      DigitizedTof->Add(digtof_p);
      digtof_p->SetSlatno(is+1);

      tup += a* fParamsTof_p->GetSigmaTime();
      time= tup/fParamsTof_p->GetTdcConv()
	     +fParamsTof_p->GetTdcOffset();
      digtof_p->SetTdcUp(time);
      
      tdown += a* fParamsTof_p->GetSigmaTime();
      time = tdown/fParamsTof_p->GetTdcConv()
	     +fParamsTof_p->GetTdcOffset();
      digtof_p->SetTdcDown(time);

	  adc_up   = 0.0;
	  adc_down = 0.0;
	 
      for(n=0;n<nlookup;n++){
	    hit_p = hittable[n];
		if(hit_p->Tof() < tup + 100){
	      ypos = hit_p->LocalYPosIn();
	      adc_up  += hit_p->Dedx()*TSonataMath::Exp(-(length-ypos)/decaylength);
		  adc_down +=hit_p->Dedx()*TSonataMath::Exp(-(length+ypos)/decaylength);
		}
      }
      adc  = adc_up * fParamsTof_p->GetADCGainUp()
	     / fParamsTof_p->GetMeVperMIP()*0.001;
      adc = adc*(1.0 + a * fParamsTof_p->GetSigmaAdc());
      digtof_p->SetAdcUp(adc);
      adc  = adc_down * fParamsTof_p->GetADCGainDown()
	     / fParamsTof_p->GetMeVperMIP()*0.001;
      adc = adc*(1.0 + a * fParamsTof_p->GetSigmaAdc());
      digtof_p->SetAdcDown(adc);



    } 
  }   
}

 void BrDigitizeTof::ListDetectorParameters() const {
  // List the current value of the digitization parameters on
  // standard out. 
  if( fParamsTof_p != NULL){
    cout <<    "The name is " << GetName() <<"n";
    fParamsTof_p->ListParameters();
  }
  else
    cout << "No parameters set for this Digitisation module" << endl;
}


 void BrDigitizeTof::Info() const{
	//
	// Standard info printout. This should give the nme of the module,
	// the author, and most recent revision. This should be implemented
	// extracting some of the info from the cvs -id etc.
	BrModule::Info();
	cout << "*n";
	cout << "*      Demonstration module written byn";
	cout << "*       Flemming Videbaekn";
	cout << "*n";
	cout << "**********************************************************n";
	cout << flush;


}