//	$Id: BrDigitizeTof.cxx,v 1.5 2002/08/30 18:44:35 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 BrTofDig. 
//      At present this code works for the following detectors
//      TOFW:  geant hits TFP1
//      TOF1:  geant hits TOF1
//      TOF2:  geant hits TOF2
//      TD1 :  geant hits D1TR
// 
//
////////////////////////////////////////////////////////////////
//

// Revisions:
//  April 4 1999, fv
//   Included BrParametersDBmanager for better detector parameter handling
//   Changed SetDetectorParamsDB to be private (scheduled for deletion)
//   Removed at the same time the delete fParamsTof_p which now is done by
//   the database manager.
//
//  July 17 2000, fv
//   Use the BrUnits::MeV to convert properly from GeantHits Dedx;
//   the units have been changed to GeV to be overall consistent.
//
#include "math.h"
#include <stdlib.h>
#include <string.h>
#include <iostream.h>

//
//  Root includes
//


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

#include "BrDigitizeTof.h"
#if !defined BRAT_BrGeometryDbManager
#include "BrGeometryDbManager.h"
#endif
#if !defined BRAT_BrParameterDbManager
# include "BrParameterDbManager.h"
#endif

#if !defined BRAT_BrDetectorParamsTof
#include "BrDetectorParamsTof.h"
#endif
#if !defined BRAT_BrMath
#include "BrMath.h"
#endif
#if !defined BRAT_BrDetectorVolume
#include "BrDetectorVolume.h"
#endif
#if !defined BRAT_BrEventNode
#include "BrEventNode.h"
#include "BrEvent.h"
#endif
#if !defined BRAT_BrDataTable
#include "BrDataTable.h"
#endif
#if !defined BRAT_BrGeantHit
#include "BrGeantHit.h"
#endif
#if !defined BRAT_BrTofDig
# include "BrTofDig.h"
#endif

ClassImp(BrDigitizeTof)

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

//_________________________________________________________________________
 BrDigitizeTof::BrDigitizeTof(const Text_t *Name, const Char_t *Title)
  :BrModule(Name,Title)
{
}

//_________________________________________________________________________
 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()
{
}

//____________________________________________________________________________
 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).

  //03-05-2000, D.Ouerdane: good place for that, Init() is called 
  //in the constructor

  //16-05-2000, Dj.O: no need of that
  //BrGeometryDbManager *gGeomDb = BrGeometryDbManager::Instance();
  //fVolumeParamsTof_p = 
  //(BrDetectorVolume*)gGeomDb->GetDetectorVolume("BrDetectorVolume", 
  //                                               (Char_t*)GetName());
  
  BrParameterDbManager *gParamDb = BrParameterDbManager::Instance();
  fParamsTof_p = (BrDetectorParamsTof*) 
    gParamDb->GetDetectorParameters("BrDetectorParamsTof",GetName());
  

}
//____________________________________________________________________________
 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.
  //
  //---------------------

  //
  // Decode the inputtable
  // The tof hits come from a table with e.g. the name "GeantHits TOF1"
  // The geant hits for tofw and Td1 is not that of the detector so the
  // selection code redefines the name.
  //
  Int_t NumHits;
  BrDataTable* GeantHits;
  BrGeantHit* ghit_p;
  BrTofDig *digtof_p;

  Char_t DigTableName[32];
  Char_t TableName[32];
  sprintf(TableName,"GeantHits %s", GetName());
  GeantHits = InputTable->GetDataTable(TableName);
  if(!GeantHits){
    if(!strcmp(GetName(),"TOFW")){
      sprintf(TableName,"GeantHits %s", "TFP1");
      GeantHits = InputTable->GetDataTable(TableName);
    }
    if(!strcmp(GetName(),"TD1")){
      sprintf(TableName,"GeantHits %s", "D1TR");
      GeantHits = InputTable->GetDataTable(TableName);
    }
  }


  if (DebugLevel() > 0)
    {
      cout << "Inside BrDigitizeTof::Event of " << GetName() << endl;
      InputTable->ListObjects();
    }
  
  if(GeantHits == 0 ) 
    {
      if(DebugLevel() > 0)
	cout << "No hits for " << GetName() << endl;
      return;
    }


  //Int_t tfp = strcmp(GetName(),"TFP1"); //for tofw panels
  NumHits = GeantHits->Entries();
  
  if(NumHits)
    {
      Int_t panel = 0;
      BrDetectorParamsTof *par;
      BrDataTable *DigitizedTof;
      sprintf(DigTableName,"DigTof %s",GetName());
//       if (!strcmp(GetName(),"TFP1"))
// 	sprintf(DigTableName,"DigTof TOFW");

      DigitizedTof = new BrDataTable(DigTableName);

      par = GetDetectorParamsTof();
      if(par == NULL)
	{
	  cout << "Parameters must be set " << endl;
	  return ;
	}

      if( DebugLevel() > 2)
	{
	  par->ListParameters();
	}
      
      // note: for each tofw panel, same set of parameters,
      // no need to worry about that
      Int_t noslats = fParamsTof_p->GetNoSlats();
      Int_t ip;     //
      Int_t slat;
      Int_t ih;     // Local hit identifier.
      
      Int_t* slattable = new Int_t [noslats];
      BrGeantHit** hittable  = new BrGeantHit* [NumHits];
      for(ip = 0; ip < noslats; ip++) slattable[ip] = 0;
      
      
      if(DebugLevel() > 0)
	cout << "Number of Input hits = " << NumHits 
	     << " in GeantHits " << GetName() << endl;
      //
      // Loop over all Slats;
      // Check on hits and create TofDig according to the algorithm given
      // in the class description.
      //

      BrGeantHit* hit_p;
      Int_t hitno = 0, nlookup, n;
      Float_t adc, 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)
		{
		  cerr << "**error ought not to happen" << endl;
		  cerr << " problem with " << GetName() << endl;
		  cerr << " number of geant hits = " << NumHits <<
		    endl;
		  cerr << " increment = " << ih << endl;
		  break;
		}
	      
	      slat = ghit_p->Isub()-1;
	      //D.Ouerdane 11-05-00: special for mrs tofw  panels
// 	      if (!tfp)
// 		{
// 		  panel = CheckSlatNo(ghit_p->Isub());
// 		  slat = (ghit_p->Isub() - panel*noslats) - 1;
// 		}
	      
	      if(DebugLevel() > 1)
		cout << " TOF name: " << GetName()
		     << "n\tSlatno: " << ghit_p->Isub()
		     << "n\tTof: " << ghit_p->Tof()
		     << "n\tTrackno: " << ghit_p->TrackNo()
		     << "n\tDedx: " << ghit_p->Dedx()
		     << "n\tYpos: " << ghit_p->LocalYPosIn()
		     << "n\tpanel slat num: " << slat+1
		     << "n\tincrement is: " << is
		     << endl;

	      if (slat == is)
		if(ghit_p->Dedx() > 0.0)
		  {
		    if(DebugLevel() > 1)
		      cout << "Add slat " << slat+1 << " i.e. Isub = "
			   << ghit_p->Isub() 
			   << " for " << GetName() << 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 BrTofDig();
	      DigitizedTof->Add(digtof_p);
	      digtof_p->SetSlatno(is+1);
	     //  if (!tfp) // for mrs tofw
// 		digtof_p->SetSlatno(is + panel*noslats + 1);
	      
	      hit_p = hittable[0];
	      
	      length = fParamsTof_p->GetScintLength();
	      dedx = hit_p->Dedx();
	      ypos = hit_p->LocalYPosIn();
	      decaylength = fParamsTof_p->GetAttLength();
	      dl   = BrMath::Exp(-(length-ypos)/decaylength);
	      // Note that Geant gives energies in GeV so instead of
	      //      adc  = dedx * fParamsTof_p->GetADCGainUp()
	      //	     / fParamsTof_p->GetMeVperMIP()*0.001*dl;
	      // we should have IGB 18FEB99
	      // Later...I should have realized that our dedx (which is
	      // actually energy deposited in the scintillator) is in MeV, and
	      // thus, the 0.001 should be removed from the above
	      // cout << "***********DEDX*******************" << dedx << endl;
	      adc  = (dedx * fParamsTof_p->GetADCGainUp() * dl)
		/ fParamsTof_p->GetMeVperMIP()/BrUnits::MeV;
	      //
	      // 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());
	      // we should have and integer adc for BrTofDig IGB 18FEB99
	      // cout << "***********adc_up #1*******************" << adc<< endl;
	      digtof_p->SetAdcUp(Int_t(adc));
	      dl   = BrMath::Exp(-(length+ypos)/decaylength);
	      adc  = dedx * fParamsTof_p->GetADCGainDown() *dl
		/ fParamsTof_p->GetMeVperMIP()/BrUnits::MeV;
	      adc = adc*(1.0 + b * fParamsTof_p->GetSigmaAdc());
	      digtof_p->SetAdcDown(Int_t(adc));
	      //cout << "***********adc_down #1*******************" << adc<< endl;
	      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(Int_t(time));
	      
	      t0  = hit_p->Tof() + (length + ypos) * 
		fParamsTof_p->GetSignalSpeed();
	      t0 += b * fParamsTof_p->GetSigmaTime();
	      time= t0/fParamsTof_p->GetTdcConv() + 
		fParamsTof_p->GetTdcOffset();
	      
	      digtof_p->SetTdcDown(Int_t(time));
	      if(DebugLevel() > 2)
		cout << "Single Hit Digitize: " 
		     << "n\tSlat no: " << hit_p->Isub() 
		     << "n\t    Tof: " << hit_p->Tof()
		     << "n\t   Dedx: " << hit_p->Dedx() 
		     << "n\t   Ypos: " << ypos << endl;
	    } 
	  //
	  //  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;
		  if(DebugLevel()>2)
		    cout << "Multihit Digitize   " <<is+1 << "   " 
			 << 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 BrTofDig 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 BrTofDig();
	      digtof_p->SetSlatno(is+1);
// 	      if (!tfp) // for mrs tofw
// 		digtof_p->SetSlatno(is + panel*noslats + 1);

	      gRandom->Rannor(a,b);

	      tup += a * fParamsTof_p->GetSigmaTime();
	      time= tup/fParamsTof_p->GetTdcConv() + 
		fParamsTof_p->GetTdcOffset();
	      digtof_p->SetTdcUp(Int_t(time));
	      
	      tdown += b* fParamsTof_p->GetSigmaTime();
	      time = tdown/fParamsTof_p->GetTdcConv() + 
		fParamsTof_p->GetTdcOffset();
	      digtof_p->SetTdcDown(Int_t(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()*BrMath::Exp(-(length-ypos)/
						  decaylength);
		      adc_down
			+=hit_p->Dedx()*BrMath::Exp(-(length+ypos)/
						    decaylength);
		    }
		}
	      // IGB 18FEB99 
	      // I have removed a factor 0.001 from adc value 
	      // see single hit stuff for explanation
	      
	      gRandom->Rannor(a,b);

	      adc  = adc_up * fParamsTof_p->GetADCGainUp() / 
		fParamsTof_p->GetMeVperMIP()/BrUnits::MeV; 
	      
	      adc = adc*(1.0 + a * fParamsTof_p->GetSigmaAdc());
	      
	      digtof_p->SetAdcUp(Int_t(adc));
	      adc  = adc_down * fParamsTof_p->GetADCGainDown() / 
		fParamsTof_p->GetMeVperMIP()/BrUnits::MeV; ;
	      adc = adc*(1.0 + b * fParamsTof_p->GetSigmaAdc());
	      
	      digtof_p->SetAdcDown(Int_t(adc));
	      DigitizedTof->Add(digtof_p);
	      
	    }
	}
      //03-05-2000, D.Ouerdane: adding table if it is not empty
      if (DebugLevel() > 0)
	{
	  cout << "tEntries of " 
	       << DigitizedTof->GetName() << ": " 
	       << DigitizedTof->GetEntries() << endl;
	}

      if (DigitizedTof->GetEntries())
	OutputTable->AddDataTable(DigitizedTof);
    }
  
  else 
    if (DebugLevel() > 0)
      cout << "tThe Geant Hit table GeantHits " 
	   << GetName() << " is empty!" << endl;
  
  if (DebugLevel() > 0)
    OutputTable->ListObjects();

}

//_________________________________________________________________________
 Int_t BrDigitizeTof::CheckSlatNo(Int_t slat) 
{

  Int_t panel = 0;
  // check slat number in tofw to give a panel id. 
  if(slat > 0 && slat <= 20)
    panel = 0;
  else if(slat > 20 && slat <= 41)
    panel = 1;
  else if(slat > 41 && slat <= 62)
    panel = 2;
  else if(slat > 62 && slat <= 83)
    panel = 3;
  else if(slat > 83 && slat <= 104)
    panel = 4;
  else if(slat > 104 && slat <= 125)
    panel = 5;

  return panel;
}

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


//_________________________________________________________________________
 void BrDigitizeTof::Print(Option_t* option) const
{
  // Standard information printout. 
  // 
  // Options: See BrModule::Print
  // 
  TString opt(option);
  opt.ToLower();
  
  BrModule::Print(option);
  if (opt.Contains("d"))
    cout << endl
	 << "  Demonstration module written by" << endl
         << "  Original author: Kris Hagel" << endl
         << "  Revisted by:     $Author: videbaek $" << endl
         << "  Revision date:   $Date: 2002/08/30 18:44:35 $"   << endl
         << "  Revision number: $Revision: 1.5 $ " << endl
         << endl
         << "*************************************************" << endl;
}

//	$Log: BrDigitizeTof.cxx,v $
//	Revision 1.5  2002/08/30 18:44:35  videbaek
//	Add capability to be able to digitize TD1 detector, which is just like a
//	the other TOF systems.
//	
//	Revision 1.4  2002/04/01 15:05:44  ejkim
//	Fixed double-use random number for TDC up and down.
//	
//	Revision 1.3  2001/08/12 13:33:02  cholm
//	Changed CTOR prototype to comply with BrModule CTOR (const Char_t*'s instead
//	of just Char_t*'s).
//	
//	Revision 1.2  2001/06/22 17:40:10  cholm
//	Changes t odo with data classes renaming.
//	
//	Revision 1.1.1.1  2001/06/21 14:55:06  hagel
//	Initial revision of brat2
//	
//	Revision 1.16  2001/03/07 12:25:27  cholm
//	* Made the method BrModule::Info() const obsolete in favour of
//	  BrModule::Print(Option_t* option="B") const. Update in other classes.
//	
//	Revision 1.15  2000/11/22 21:17:50  videbaek
//	Moved parameter and geometry setup to Init() and/or removed Init() from constructor.
//	
//	Revision 1.14  2000/09/13 16:41:28  ouerdane
//	Removed TFPx panel tables. TOFW in simulation is now treated like raw data without any explicit splitting
//	
//	Revision 1.13  2000/09/04 15:03:20  videbaek
//	Remove include reference to <unistd.h> that does not exist for WIN32
//	and is not needed.
//	
//	Revision 1.12  2000/07/17 19:51:08  videbaek
//	Change to GeV as unit for momentum and energy.
//	
//	Revision 1.11  2000/05/17 10:28:50  ouerdane
//	see top of file
//	
//	Revision 1.11  2000/05/14 ouerdane
//	Now can digitize mid rapidity TOFW panels.
//	
//	Revision 1.10  2000/04/06 20:27:00  cholm
//	Corrected some minor faults
//	
//	Revision 1.9  1999/07/03 21:06:49  hagel
//	Change output table name
//	
//	Revision 1.8  1999/04/14 20:10:40  videbaek
//	Utilize DbManager to get defailt parameters
//
//	Revision 1.7  1999/02/19 21:32:40  bearden
//	Added ADCs to Tof calibration, now returns ADC sum
//	(top+bottom) in channels
//	Changed the way the hodoscope hit position is compared to the
//	track position.  This is now done in the "local coordinates" 
//      of the hodo.
//
//      
//	Revision 1.6  1999/01/28 21:26:13  videbaek
//	Inserted DebugLevel protection around some debug output
//
//	Revision 1.5  1999/01/28 16:33:15  hagel
//	Remove references to TSonataMath
//
//	Revision 1.4  1998/12/21 20:28:47  videbaek
//	updated makeNT file. Cosmetics for tof
//
//	Revision 1.3  1998/09/08 15:31:44  alv
//	Added fVolumeParamsTof_p=0 in constructor.
//
//	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
//

Copyright ; 2002 BRAHMS Collaboration <brahmlib@rcf.rhic.bnl.gov>
Last Update on by

Validate HTML
Validate CSS