//   $Id: BrDigitizeTPC.cxx,v 1.10 1998/08/05 19:34:27 hagel Exp $
//      File BrDigitizeTPC.cxx
//
//      $Log: BrDigitizeTPC.cxx,v $
//      Revision 1.10  1998/08/05 19:34:27  hagel
//      Fix BrDigTPC memory leak
//
//      Revision 1.9  1998/07/02 13:58:43  videbaek
//       comment updates mainly
//
//      Revision 1.8  1998/05/14 21:18:15  videbaek
//      update makefile - do not reload .so libraries
//
//      Revision 1.7  1998/05/13 20:30:50  hagel
//      Changes to compile under Solaris
//
//      Revision 1.6  1998/05/08 20:45:56  videbaek
//      TPC class updates; added functionality
//
//      Revision 1.5  1998/05/08 00:58:43  videbaek
//      TPC code updates
//
//      Revision 1.4  1998/04/30 16:40:48  hagel
//      Changes to work with new version of GBRAHMS
//
//      Revision 1.3  1998/04/06 21:12:00  videbaek
//      Clean up and additions for Win95
//
//	Revision 1.2  1998/04/01 22:08:36  hagel
//	Add Local tracking and fix bugs in Digitization
//
//	Revision 1.1.1.1  1998/03/04 21:34:03  brahmlib
//	Brat tpc
//
//

////////////////////////////////////////////////////////////////
//
// The BRDigitize TPC is the digitization module for the
// pad TPC's. The module performs a slow digitization using the
// geant hit information.
// 
////////////////////////////////////////////////////////////////

//
//	System includes
//
#include "math.h"
#include <stdlib.h>
#ifdef UNIX
#include <unistd.h>
#endif
#include <iostream.h>
//
//  Root includes
//
#include "BrDigitizeTPC.h"
#include "btwn.h"
#include "TMath.h"
#include "TSonataMath.h"

#ifndef ROOT_TRandom
#include "TRandom.h"
#endif
#include "TH1.h"


#include "BrDetectorParamsTPC.h"
#include "BrDetectorVolume.h"
#include "BrEventNode.h"
#include "BrDigitizeTPC.h"
#include "BrDataTable.h"

#include "BrGeantHit.h"
#include "BrDigTPC.h"

ClassImp(BrDigitizeTPC)

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

  BrDigitizeTPC::BrDigitizeTPC(Text_t *Name,Char_t *Title)
   : BrModule(Name, Title)
{
//padres_p = new BrClonesArray("BrDigTPC",10000);
for(Int_t i=0;i<max_irow;i++) padres[i] = new BrDataTable();
}

 void BrDigitizeTPC::SetDetectorParamsTPC(BrDetectorParamsTPC* par)
{ 
 // Set the TPC parameters with the values of the external parameters
 // referenced by pointer
 // Note that no owner ship is transferred. A copy is made.
 // 
 // 
	
	if( fParams_p != 0)
		delete fParams_p;
	fParams_p = new BrDetectorParamsTPC("Dummy","TPCINT");
	*(BrDetectorParamsTPC*)fParams_p = *par;
	// The conversion cast is needed because fParams_p is of
	// type BrDetectorParams.
}

 void BrDigitizeTPC::SetDetectorParamsTPC(const BrDetectorParamsTPC& par)
{
  // Set the TPC 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( fParams_p != 0)
		delete (BrDetectorParamsTPC*) fParams_p;
	fParams_p = new BrDetectorParamsTPC("Dummy","ParamsTPC");
	*(BrDetectorParamsTPC*)fParams_p = par;
	// The conversion cast is needed because fParams_p is of
	// type BrDetectorParams.
}

 void BrDigitizeTPC::ListDetectorParameters(){
// List the current value of the digitization parameters on
// standard out. 
if( fParams_p != NULL)
	fParams_p->ListParameters();
else
cout << "No parameters set for this Digitisation module" << endl;
}


 void BrDigitizeTPC::EventStatisticsStart() {
  fTimer.Start();
}

 void BrDigitizeTPC::EventStatisticsEnd() {
  fTimer.Stop();
  fCpuTime += fTimer.CpuTime();
}

 void BrDigitizeTPC::ListEventStatistics() {
  cout<<"Accumulated CPU time in "<<GetName()<<" is "<<fCpuTime<<endl;
  cout<<"NumGeantHits = "<<fNumGeantHits<<" NumDigTPC = "<<fNumDigTPC<<endl;
}

 BrDigitizeTPC::~BrDigitizeTPC(){
  //
  // Is there a problem with the fact the the fParams_p can point to different kinds
  // of datastructures i.e. potential memory leak. probably yes since fParams_p is
  // defined as BrDetectorParams nor BrDetectorParamsTPC.
  // By second thought it is probably ok. As long as all parameters are defined/set
  // Using the SetDetectorParamsTPC / DC /BB etc..
  //
  if( fParams_p != 0)
    delete fParams_p;
  if( fVolumeParams_p != 0)
    delete fVolumeParams_p;
  
}

 void BrDigitizeTPC::Init(){
  // The init member function/method is used to define histograms
  // The histograms is added to the list of hist objects in this
  // module.
  d_Hhits = new TH1F("Hhits","Geant Hits per event",40, (float)0., (float)40.);
  HistogramList()->Add(d_Hhits);
}

 void BrDigitizeTPC::Event(BrEventNode* InputTable, BrEventNode* OutputTable)
{
  // Event method to be called once per event. This performes the actual slow simulation
  // of the hits on the TPC.
  //
  Bool_t hyperbolic = kFALSE;
  
  const double pi = 3.14159265358979323846;
  const Float_t  twopi = (Float_t)2.0 * (Float_t)pi;
  const Float_t  sqpi = (Float_t)sqrt(pi);
  const Float_t  sqtwopi = (Float_t)sqrt(2 * pi);
  const Float_t  sqtwo = Float_t(sqrt(2.));
  const Float_t  dedx_to_n0 = (Float_t)(1.0 / .018);
  
  //following changed from 1000 10-7-97 kh
#define max_padres 5000		//50000 in fortran version
//#define max_irow   20
  
/*  struct {
    Int_t padno;
    Int_t time;
    BrGeantHit *ghits;
    Int_t used;
    Float_t adc;
  } padres[max_padres][max_irow];

  Int_t npadres[max_irow];*/
  Int_t npadres;
//  BrDataTable padres[max_irow];
  
  Float_t a,b,z,ax,ay,n0,xlo,xhi,dxl,dxt,xc,yc;
  Float_t xbase,ybase,time_bucket,Anode_Gap,sigma_0,sigma_prf,xb,path_length;
  Float_t cm_to_timeb,dedx_cm,factor_to_channels,xwire,ywire;
  Float_t sigma_T,sigma_L,diffuse_T,diffuse_L,Ldrift;
  Float_t y_low,y_high,tnorm,pnorm,edge_low,edge_high,qpad,qbin;
  
  Int_t i,ip;
  Int_t irow,ihit,ixmin,ixmax,min_time,max_time,ipad,itime;
  
  Float_t detZ;
  
  BrDetectorParamsTPC *par;
  BrDetectorVolume *vol;
  Int_t NumHits;
  BrDataTable* GeantHits;
  BrGeantHit* ghit_p;
  BrDigTPC *digtpc_p;
  BrDigTPC *digtpc_p1;
  BrDigTPC *digtpc_p2;

  if(DebugLevel()>0){ 
    printf("Inside Digitize of %sn",GetName());;
    InputTable->ListObjects();
  }
  vol = GetDetectorVolume();
  if(vol == NULL){
    cout << GetName() << "Volume must be set " << endl;
    return ;
  }
  if(DebugLevel()>0 ) vol->ListParameters();
  par = GetDetectorParamsTPC();
  if(par == NULL){
    cout << "Parameters must be set " << endl;
    return ;
  }
  if( DebugLevel() > 2){
    cout << " Getting Parameters" << endl;
    par->ListParameters();
  }
  // Decode the inputtable
  //
  Char_t TableName[32];
#ifndef OLD_CDAT_VERSION
  sprintf(TableName,"GeantHits %s\x0", GetName());
#else
  sprintf(TableName,"GeantHits %sPR\x0", GetName());
#endif
  GeantHits = InputTable->GetDataTable(TableName);
  
  if(GeantHits == 0 ) {
    if(DebugLevel() > 0){
      cout << "No Table " << GetName() << endl;
    }
    return;
  }
  
  NumHits = GeantHits->GetEntries();
  if(NumHits == 0 ) {
    if(DebugLevel() > 0){
      cout << "No Hits in " << GetName() << endl;
    }
    return;
  }
  //
  // ObjectList returns TObjectArray
  // Prepare the output table
  //
  sprintf(TableName,"DigitizedTPC %sx0",GetName());
  BrDataTable *DigitizedTPC = new BrDataTable(TableName);
  OutputTable->AddObject(DigitizedTPC);
  
//  if( NumHits > 50 ) return;	//too many hits for comfort
  
  
  //detZ = DetectorElementPosition[idet];
  detZ = (Float_t)0.0;	//For now with TPC for test
  
  if( vol->GetCut() == 0 ) {
    ybase = -vol->GetSize()[1] / (Float_t)2.;
    xbase = -vol->GetSize()[0] / (Float_t)2.;
  }
  else {
    ybase = -vol->GetFsize()[1] / (Float_t)2.;
    xbase = -vol->GetFsize()[0] / (Float_t)2.;
  }
  
  time_bucket = par->GetTimeBucket();
  Anode_Gap   = par->GetAnodeGap();
  sigma_0     = par->GetShapingTime() /2.36;
  sigma_prf   = par->GetAnodeGap();

  if(DebugLevel() > 1) {
    cout << " time  bin " << time_bucket << endl;
    cout << " Anode gap " << Anode_Gap << endl;
    cout << " Sig(PRF)  " << sigma_prf << endl;
    cout << " Sig(time)  " << sigma_0 << endl;
  }
  xb = -par->GetPadsprow() * par->GetPadDistance();
  
  xbase = TMath::Max(xb,xbase);
  
  cm_to_timeb = (Float_t)1.0/(par->GetDriftv() * time_bucket);
  factor_to_channels = par->GetADCGain() * par->GetADCChannels() / (Float_t)100.;

  
  for(irow=0;irow<max_irow;irow++) padres[irow]->Clear();

//for(irow=0;irow<par->GetNumberOfRows();irow++) {
//   npadres = 0;
    for(ihit=0;ihit<NumHits;ihit++) {
      ghit_p = (BrGeantHit*)GeantHits->At(ihit);
	  irow = ghit_p->Isub() - 1;		//-1 because padrow number starts with 1
      ax = ( ghit_p->LocalPosOut()[0] - ghit_p->LocalPosIn()[0] ) / 
	( ghit_p->LocalPosOut()[2] - ghit_p->LocalPosIn()[2] );
      ay = ( ghit_p->LocalPosOut()[1] - ghit_p->LocalPosIn()[1] ) / 
	( ghit_p->LocalPosOut()[2] - ghit_p->LocalPosIn()[2] );


//    The calculation of path_length probabally has to be revisited.  In the
//    previous version the positions were the front and back of the TPC.  Now it
//    is the front and back of the padrow.  On the other hand, the dedx also seems
//    to be scaled.  dedx_cm gives numbers of the same order as the ones in the 
//    previous version.
      path_length = (Float_t)0.0;
      for(i=0;i<3;i++)path_length+=(Float_t)pow(ghit_p->LocalPosOut()[i]-ghit_p->LocalPosIn()[i],2.);
      path_length = (Float_t)sqrt(path_length);
      dedx_cm = ghit_p->Dedx() / path_length;
      n0 = dedx_cm * dedx_to_n0;
//    z = detZ + par->GetRowDistance() * ( irow+1-(par->GetNumberOfRows()+1)/(Float_t)2. );
//    I am assuming that LocalPosIn and Out is the position of the submodule
//    Get the z of the center of the padrow.
      z = (ghit_p->LocalPosOut()[2] + ghit_p->LocalPosIn()[2]) / (Float_t)2.0;
      dxt = par->GetPadLength() * TMath::Abs(ax) / (Float_t)2.;
      dxl = par->GetPadLength() * TMath::Abs(ay) / (Float_t)2.;
#ifdef debug
	  cout << " position" << z << " " << dxt << " " << dxl << endl;
#endif

      if( BTWN(ghit_p->LocalPosIn()[2],z,ghit_p->LocalPosOut()[2]) ) {
         xc = ghit_p->LocalPosIn()[0] + ax * ( z - ghit_p->LocalPosIn()[2] );
         yc = ghit_p->LocalPosIn()[1] + ay * ( z - ghit_p->LocalPosIn()[2] );

         Ldrift = yc - ybase;
         diffuse_T = par->GetDtrans() * (Float_t)TMath::Sqrt( Ldrift+2. );
         diffuse_L = par->GetDlong()  * (Float_t)TMath::Sqrt( Ldrift+2. );

/*       Use a 3 sigma limit on searching for pads which might have been hit
         At this point we have all the cluster variables needed
         ax = tan(alpha), ay=tan(lambda)
         yc = Ldrift, xc which are positions before evaluating the diffusion
         in drifts and response widths
*/
         n0 *= par->GetPadLength() * (Float_t)sqrt( 1.0 + ax*ax + ay*ay );

//       Evaluate the position on the Anode wires..

         xwire = xc - xbase;
         ywire = Ldrift;

//       Distribute signal according to n0 on the Pad

         gRandom->Rannor(a,b);
         xwire += diffuse_T / (Float_t)TMath::Sqrt(n0) * a;
         ywire += diffuse_L / (Float_t)TMath::Sqrt(n0) * b;

//       Evaluate the pad response
//
         xlo = xwire - (Float_t)3.0 * Anode_Gap;
         xhi = xwire + (Float_t)3.0 * Anode_Gap;

	 if( DebugLevel() > 2){
	   cout.width(10);
	   cout << "     DT " << diffuse_T << " DL" << diffuse_L << endl;
	   cout << "     n0 " << n0 << endl;
	   cout << "     " << xlo << " " << xhi << endl;
	 }

//       Find pads which must be looked at.

         ixmin = (Int_t)( xlo / par->GetPadDistance() ) + 1;
         ixmin = TMath::Max(ixmin,1);
         ixmax = (Int_t)( xhi / par->GetPadDistance() ) + 1;
         ixmax = TMath::Min(par->GetPadsprow(),ixmax);

/*       Pad numbers are counted 1,2,3.. from lower X values. 
         there is the same number of Pads on each side of the centerline

         longitudinal Diffusion i.e. along Drift
         Convert into time buckets

         Do fold the response of the shaper into this
*/
         sigma_L = (Float_t)sqrt(pow(sigma_0*par->GetDriftv(),2) + diffuse_L*diffuse_L);
         sigma_T = (Float_t)sqrt(sigma_prf*sigma_prf             + diffuse_T*diffuse_T);

/*       As a temporary solution to the problem of angled tracks, fold in the
         sigma of the square response of ax and ay
*/
         sigma_L = (Float_t)sqrt( sigma_L*sigma_L + dxl*dxl/(Float_t)12.);
         sigma_T = (Float_t)sqrt( sigma_T*sigma_T + dxt*dxt/(Float_t)12.);

         y_low  = ywire - (Float_t)3.0 * sigma_L;
         y_high = ywire + (Float_t)3.0 * sigma_L;
         min_time = (Int_t)(y_low * cm_to_timeb);
         min_time = TMath::Max(min_time,1);
         max_time = (Int_t)(y_high * cm_to_timeb);
         max_time = TMath::Min(max_time,256);
         tnorm = (Float_t)1. / ( sqtwopi * sigma_L );
         pnorm = (Float_t)1. / ( sqtwopi * sigma_T );
	 if( DebugLevel() > 2){
	   cout << "Padrow #" << irow <<endl;
	   cout.width(10);
	   cout << "     " << ixmin << " "<< ixmax << endl;
	 }
	 if( DebugLevel() > 3){
	   cout.width(10);
	   cout << "     " << min_time << " " << max_time << endl;
	 }
	 
         for( ipad = ixmin;ipad<=ixmax;ipad++ ) {
	   edge_low  = ((Float_t)ipad - (Float_t).5)*par->GetPadDistance() - par->GetPadWidth()/(Float_t)2.;
	   edge_high = ((Float_t)ipad - (Float_t).5)*par->GetPadDistance() + par->GetPadWidth()/(Float_t)2.;
	   if( hyperbolic ) {
	     printf("This part has not been tested yet!!!");
	     qpad = (Float_t)atan((Float_t)sinh(pi*(edge_high-xwire)/(Float_t)2./Anode_Gap))
	       - (Float_t)atan((Float_t)sinh(pi*(edge_low -xwire)/(Float_t)2./Anode_Gap));
	   }
	   else {
	     qpad = pnorm*(Float_t)0.5*(TSonataMath::Erf((edge_high-xwire)/(sigma_T*sqtwo))
					- TSonataMath::Erf((edge_low -xwire)/(sigma_T*sqtwo)));
	   }
	   
	   for( itime = min_time;itime<=max_time;itime++ ) {
	     edge_low  = (Float_t)itime     / cm_to_timeb;
	     edge_high = (Float_t)(itime+1) / cm_to_timeb;
	     qbin = (Float_t)0.5 * (TSonataMath::Erf((edge_high - ywire) / (sigma_L*sqtwo))
				    -  TSonataMath::Erf((edge_low  - ywire) / (sigma_L*sqtwo)));
	     
	     qbin *= tnorm * qpad * n0;
	     
	     if( qbin > .05 ) {
	       digtpc_p = new BrDigTPC();
	       padres[irow]->Add(digtpc_p);
	       
	       digtpc_p->SetIrow(irow);
	       digtpc_p->SetPadno(ipad);
	       digtpc_p->SetTime(itime);
	       digtpc_p->SetAdc(qbin * factor_to_channels);
	       digtpc_p->SetUsed(0);
	       digtpc_p->SetNhits(1);
/*			      if(npadres<max_padres) {
                     padres[npadres].padno = ipad;
                     padres[npadres].time  = itime;
                     padres[npadres].adc   = qbin * factor_to_channels;
                     padres[npadres].ghits = ghit_p;
                     padres[npadres].used  = 0;
                     npadres++;
				     }
				  else {
				     printf("npadres is getting too large, irow is %dn",irow);
				     }*/
	     }	//qbin test
	   }	//loop over itime
	 }		//loop over ipad
      }			//if BTWN
    }				//loop over ihit
/* 
   We have now generated the 'raw' hits.  They should now be merged for
   overlapping hits. One digTPC will be created for each hit found.
   We need to be sure that these all  deleted when the event is
   done or else we will have a severe memory leak.
*/


    for(irow=0;irow<par->GetNumberOfRows();irow++) {
      
      npadres = padres[irow]->GetEntries();
      if( DebugLevel() > 1){
	cout << "Padrow #" << npadres << " entries" <<endl;
      }
  
      for(ip=0;ip<npadres;ip++) {
	digtpc_p1 = (BrDigTPC*)padres[irow]->At(ip);
	if( digtpc_p1->GetUsed() == 0 ) {
	  
	  digtpc_p = new BrDigTPC();
	  DigitizedTPC->Add(digtpc_p);
	  
	  digtpc_p->SetIrow(irow);
	  digtpc_p->SetPadno(digtpc_p1->GetPadno());
	  digtpc_p->SetTime(digtpc_p1->GetTime());
	  digtpc_p->SetAdc(digtpc_p1->GetAdc());
	  digtpc_p->SetUsed(0);
	  digtpc_p->SetNhits(1);
	  
	  for(Int_t ip2=ip+1;ip2<npadres;ip2++) {
	    digtpc_p2 = (BrDigTPC*)padres[irow]->At(ip2);
	    if( digtpc_p1->GetPadno() == digtpc_p2->GetPadno() &&
		digtpc_p1->GetTime()  == digtpc_p2->GetTime() ) {
	      //	        if( padres[ip][irow].padno == padres[ip2][irow].padno &&
	      //	           padres[ip][irow].time  == padres[ip2][irow].time ) {
	      //			   digtpc_p->AddAdc(padres[ip2][irow].adc);
	      digtpc_p->AddAdc(digtpc_p2->GetAdc());
	      
	      //	   Need to put in histogram generation here
	      //	           padres[ip2][irow].used = 1;
	      digtpc_p2->SetUsed(1);
	      digtpc_p->IncNhits();
	      
	      //     NEED TO  THIS TO USE ROOT TABLES???
	      //     IT IS NOT USED IN SONATA OTHER THAN TO BE CREATED!!!
	      //     This are used in the old analysis to check e.g. purity of found
	      //	   tracks checking how real to fiund correlates. See the old analysis.F
	      //	   fv comment
	      //	   tpchit[NumTPCHit] = new tpcHIT();
	      //	   tpchit_p = tpchit[NumTPCHit];
	      //	   NumTPCHit++;
	      //	   if(NumTPCHit > maxtpchit) cout<<"NumTPCHit too large"<<endl;
	      
	      //	   tpchit_p->Imod = irow + 1;  //+1 is to match fortran
	      //	   tpchit_p->digtpc = digtpc_p;
	      //	   tpchit_p->ghits  = padres[ip2].ghits;
	      
	    }	        //if padres...
	  }		//loop over ip2
	}               //if padres...
      }                 //loop over ip
    }			//loop over irow
    
//  for(irow=0;irow<par->GetNumberOfRows();irow++) 
    for(irow=0;irow<max_irow;irow++) 
      padres[irow]->Delete();
    
//Get stuff for statistics

    fNumGeantHits = GeantHits->GetEntries();
    fNumDigTPC = DigitizedTPC->GetEntries();

}


/*
  void BrDigitizeTPC::LocalTrack(BrEvent* InputTable,BrEvent* OutputTable)
  {
  Int_t status;
  
  //These next three need to be loaded in.  Are they done in LoadCuts??
  static Float_t dx = (Float_t)1.5;
  static Float_t dy = (Float_t)0.6;
  static Int_t   maxmiss = 2;
  //LoadCuts();????? Appears never to come there in FORTRAN Sonata
  
  CombinedHits->Clear();
  LocalTracks->Clear();
  TrackHits->Clear();
  VirtualTracks->Clear();
  
  TPCPadCluster(InputTable);
  status = TPCPadTrack(0,maxmiss,dx,dy);
  if(!status) FillLocalTracks(OutputTable);
*/