//  $Id: BrDigitizeDC.cxx,v 1.12 1998/07/03 15:47:29 hagel Exp $
//
// $Log: BrDigitizeDC.cxx,v $
// Revision 1.12  1998/07/03 15:47:29  hagel
// Clean up g++ warnings
//
// Revision 1.11  1998/07/02 17:54:15  hagel
// Add const to GetDetectorVolume
//
// Revision 1.10  1998/07/01 21:26:03  videbaek
//  fixed some comments
//
//
#include "BrDigitizeDC.h"
#include "BrDCPlane.h"
#include "BrDetectorParamsDC.h"
#include "BrDetectorVolume.h"
#include "BrGeantHit.h"
#include "BrDigHit.h"
#include "BrEventNode.h"
#include "btwn.h"
#include "TSonataMath.h"
#include "TRandom.h"
#include "BrDCHitPosition.h"

#include <iostream.h>


ClassImp(BrDigitizeDC)

////////////////////////////////////////////////////////////////
//
// The TBrDigitizeDC is the base module for digitization of drift
//  Chambers. The algorithms are based on those in Sonata and Sonata++ 
//
//
////////////////////////////////////////////////////////////////

 BrDigitizeDC::BrDigitizeDC()
:BrModule()
{
}

 BrDigitizeDC::BrDigitizeDC(Char_t *Name,Char_t *Title)
:BrModule(Name,Title)
{
fParams_p       = NULL;
for(Int_t i=0;i<4;i++) fVolumeParams_p[i] = NULL;
fCpuTime = 0;
Init();
}


 void BrDigitizeDC::SetDetectorParamsDC(BrDetectorParamsDC* par)
{ 
 // Set the DC 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 BrDetectorParamsDC("Dummy","DCINT");
	*fParams_p = *par;
	// The conversion cast is needed because fParams_p is of
	// type BrDetectorParams.
}

 void BrDigitizeDC::SetDetectorParamsDC(const BrDetectorParamsDC& 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  fParams_p;
	fParams_p = new BrDetectorParamsDC("Dummy","ParamsTPC");
	*fParams_p = par;
	// The conversion cast is needed because fParams_p is of
	// type BrDetectorParams.
}

 void BrDigitizeDC::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 BrDigitizeDC::SetDetectorVolume(BrDetectorVolume *vol) {
Int_t i;
  for(i=0;i<fNumVol;i++) {
    if(!strcmp(fVolumeParams_p[i]->GetName(),vol->GetName())) {
       fVolumeParams_p[i] = vol;
       return;
    }
  }
if(fNumVol<4) {
   fVolumeParams_p[fNumVol] = vol;
   fNumVol++;
   }
else cout<<"Too many volume parameters, this one not added"<<endl;
}

 BrDetectorVolume *BrDigitizeDC::GetDetectorVolume(const Char_t *name) {
Int_t i;
for(i=0;i<fNumVol;i++) {
   if(!strcmp(fVolumeParams_p[i]->GetName(),name)) {
      return fVolumeParams_p[i];
      }
   }
cout<<"Detector Volume "<<name<<" does not exist"<<endl;
return 0;
}

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

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

 void BrDigitizeDC::ListEventStatistics() {
  Char_t buf[32];
  cout<<"Accumulated CPU time in "<<GetName()<<" is "<<fCpuTime<<endl;
  for(Int_t i=0;i<3;i++) {
     sprintf(buf,"NumGeantHits for %s A%d is x0",GetName(),i);
     cout<<buf<<fNumGeantHits[i]<<endl;
     }
  cout<<" NumDigHits = "<<fNumDigHits<<endl;
}

 BrDigitizeDC::~BrDigitizeDC(){
  //
  // 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;
  for(Int_t i=0;i<4;i++) {
     if( fVolumeParams_p[i] != 0) delete fVolumeParams_p[i];
     }
  
}

 void BrDigitizeDC::Init(){
  // The init member function/method is used to define histograms
  // The histograms is added to the list of hist objects in this
  // module.
Char_t HistName[32];
sprintf(HistName,"Hhits%s",GetName());
d_Hhits = new TH1F(HistName,"Geant Hits per event",40, (float)0., (float)40.);
sprintf(HistName,"DZ%s",GetName());
d_Dz    = new TH1F(HistName,"Distribution of Z",310,-15.0,15.0);
sprintf(HistName,"%sxz",GetName());
//d_Graphxz = new TGraph();

//  HistogramList()->Add(d_Hhits);
}

 void BrDigitizeDC::Event(BrEventNode* InputTable,BrEventNode* OutputTable)
  // Event method for DC digitization
  //
{
Digitize(InputTable,OutputTable);
Decode(OutputTable);
}

 void BrDigitizeDC::Digitize(BrEventNode* InputTable,BrEventNode* OutputTable)
{
static const Int_t digmode = 1;
static const Double_t rad = .01745329;
//Float_t FGauss();
Int_t i;
Int_t ihit,iplane;
Float_t ax,ay,z,zsp_eff,cosT,sinT;
Float_t r,xc,yc,uc,au,uin,uout;
Float_t xwire,u0,z0,u1,z1,u2,z2,dU;
Int_t iuin,iuout,Istep,uwire;

//TDriftChamberPlane* dcpl;
BrDCPlane* dcpl;
Int_t numplanes;

Int_t NumHits;
Int_t idet;

Float_t detZ;
Float_t SubModulePosition[3] = {(Float_t)-30.0,(Float_t)0.0,(Float_t)30.0};


BrDetectorParamsDC *par;
BrDetectorVolume *vol[4];

BrDataTable* GeantHits[3];
BrGeantHit* ghit_p;
BrDigHit *dighit_p;

//InputTable->ListObjects();
//printf("Inside Digitize of %sn",GetName());;

//vol->ListParameters();
//fParams->ListParameters();

//Get stuff from the input table
Char_t TableName[32];
for(idet=0;idet<3;idet++ ) {
   sprintf(TableName,"GeantHits %sA%dx0",GetName(),idet+1);
   GeantHits[idet] = InputTable->GetDataTable(TableName);
   if(GeantHits[idet] == 0 ) {
      if(DebugLevel() > 0){
         cout << "No hits for " << GetName() << "Module number "<<idet<<endl;
         }
      return;
      }
   }
//NumHits = GeantHits->GetEntries();
  
// GetObjectList returns TObjectArray

//Prepare the output table
sprintf(TableName,"%s DigitizedDCx0",GetName());
BrDataTable *DigitizedHits = new BrDataTable(TableName);
OutputTable->AddObject(DigitizedHits);

//if( NumHits > 50 ) { 	//too many hits for comfort
//   printf("DC: %s Rejecting this event, too many hitsn",GetName());
//   SetNoFurtherExecute();
//   return;
//   }

par = GetDetectorParamsDC();
if(par == NULL){
   cout << GetName() << "::Digitize; Parameters must be set " << endl;
   return ;
   }

vol[3] = GetDetectorVolume(GetName());
if(vol[3] == NULL){
   cout << GetName() << "::Digitize; Volume must be set " << endl;
   return ;
   }

for(idet=0;idet<3;idet++) {

NumHits = GeantHits[idet]->GetEntries();

if(DebugLevel()>0) {
   cout<<"Starting Digitization for "<<GetName()<<endl;
   for(ihit=0;ihit<NumHits;ihit++) {
      cout<<(BrGeantHit*)GeantHits[idet]->At(ihit);
      }
   }

//Initialize things to make g++ happy (or at least not sad)
au    = (Float_t)0.0;
xwire = (Float_t)0.0;
u2    = (Float_t)0.0;
z2    = (Float_t)0.0;

detZ = SubModulePosition[idet];	//Need to change!!!

d_Hhits->Fill(NumHits);

numplanes = par->GetNplane();
for( ihit=0;ihit<NumHits;ihit++) {
   ghit_p = (BrGeantHit*)GeantHits[idet]->At(ihit);
   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] );

   d_Dz->Fill(ghit_p->LocalPosIn()[2]);
   d_Dz->Fill(ghit_p->LocalPosOut()[2]);

   zsp_eff = par->GetPldis() *(Float_t)(0.8 / 2.);
   for(iplane=0;iplane<numplanes;iplane++) {
      dcpl = par->GetPlaneAt(iplane);
      cosT = (Float_t)TMath::Cos(dcpl->GetWirang() * rad);
      sinT = (Float_t)TMath::Sin(dcpl->GetWirang() * rad);
      //      z = detZ + par->GetPldis()*((Float_t)(iplane+1)-(Float_t)(numplanes+1)/(Float_t)2.);
      z = par->GetPldis()*((Float_t)(iplane+1)-(Float_t)(numplanes+1)/(Float_t)2.);
//    if(BTWN(ghit_p->LocalPosIn()[2],z,ghit_p->LocalPosOut()[2])) {
//    Following changed because LocalPosIn and Out are now relative to the submodule
//    so to see if they are in the detector, one needs to take into account position of
//    module.
      if(DebugLevel() > 4) {
	     cout<<"Digitize z position is "<<z<<endl;
         }
//    if(BTWN(detZ+ghit_p->LocalPosIn()[2],z,detZ+ghit_p->LocalPosOut()[2])) {
      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]);
	 uc = xc * cosT - yc * sinT;
	 au = ax * cosT + au * sinT;
	 uin  = uc - au * zsp_eff;
	 uout = uc + au * zsp_eff;
	 if( digmode == 1 ) {
	    uin  = uc;
	    uout = uc;
	    }
	 if( dcpl->GetNwire() == 1 ) {
//	    TMath::Int appears to behave as FORTRAN NINT
	    iuin  = TSonataMath::Int(uin /dcpl->GetWirdis());
	    iuout = TSonataMath::Int(uout/dcpl->GetWirdis());
	    }
	 else {
//	    (Int_t) appears to behave as FORTRAN INT
	    iuin  = (Int_t)(uin /dcpl->GetWirdis());
	    iuout = (Int_t)(uout/dcpl->GetWirdis());
	    if( uin < 0 ) iuin = iuin - 1;
	    if( uout< 0 ) iuout = iuout - 1;
	    }
         if(iuin == iuout) {
	    r = gRandom->Rndm();
	    if(r<= par->GetEff()) {
	       if( dcpl->GetNwire() == 1 ) {
	          xwire = (Float_t)iuin * dcpl->GetWirdis();
		  }
	       else xwire = ((Float_t)iuin + (Float_t)0.5) * dcpl->GetWirdis();

	       dighit_p = new BrDigHit();
	       DigitizedHits->Add(dighit_p);
	       dighit_p->SetID(DigitizedHits->GetEntries());
	       dighit_p->SetIdet(idet);	//+1??to be like FORTRAN
	       dighit_p->SetImod(iplane); //should it be +1 to be like FORTRAN
	       dighit_p->SetWire(iuin);
	       dighit_p->SetTime((Int_t)((TMath::Abs(uc - xwire) + par->GetPosres()*TSonataMath::FGauss())/par->GetDriftv()));
	       dighit_p->SetWidth((Int_t)(par->GetTwoPar() / par->GetDriftv()));
	       dighit_p->SetIfuse(0);
	       dighit_p->SetUsed(0);
	       dighit_p->SetDcpl(iplane);
	       }
	    }
         else {
	    cout<<"Digitize: This part not debugged, use with caution (debug before using)"<<endl;
	    if( au > 0 ) Istep = 1; else Istep = -1;
	    z0 = au*(xwire-uin+au*(-zsp_eff))/(au*au+1);
	    u0 = -z0 / au + xwire;
//	    (z0,u0): the nearest point on the path from the wire
	    for(uwire=iuin;uwire<=iuout;uwire+=Istep) {
	    if(dcpl->GetNwire()==1) xwire = (Float_t)uwire             * dcpl->GetWirdis();
	    else                    xwire =((Float_t)uwire+(Float_t)0.5) * dcpl->GetWirdis();
	    if( uwire == iuin ) {
	       u1 = uin;
	       z1 = -zsp_eff;
	       if(dcpl->GetNwire()==1) u2 = ((Float_t)iuin +                (Float_t)Istep*(Float_t)0.5)*dcpl->GetWirdis();
	       else		       u2 = ((Float_t)iuin + (Float_t)0.5 + (Float_t)Istep*(Float_t)0.5)*dcpl->GetWirdis();
	       z2 = z1 + (u2 - u1) / au;
	       }
	    else if( uwire == iuout ) {
	       u1 = u2;
	       z1 = z2;
	       u2 = uout;
	       z2 = zsp_eff;
	       }
            else {
	       u1 = u2;
	       z1 = z2;
	       u2 = u1 + (Float_t)Istep * dcpl->GetWirdis();
	       z2 = z1 + (u2 - u1) / au;
	       }
	    if(z0<-zsp_eff) {
	       dU = (Float_t)TMath::Sqrt(zsp_eff*zsp_eff+(uin-xwire)*(uin-xwire))+par->GetPosres()*TSonataMath::FGauss();
	       }
	    else if(z0>zsp_eff) {
	       dU = (Float_t)TMath::Sqrt(zsp_eff*zsp_eff+(uout-xwire)*(uin-xwire))+par->GetPosres()*TSonataMath::FGauss();
				  }
	    else if(u0<xwire-0.5*dcpl->GetWirdis()||u0>xwire+0.5*dcpl->GetWirdis()) {
	       dU = (Float_t)0.5 * dcpl->GetWirdis() + par->GetPosres() * TSonataMath::FGauss();
	       }
	    else {
	       dU = (Float_t)TMath::Abs(-xwire-au*z1+u1)/(Float_t)TMath::Sqrt(au*au+1.)+par->GetPosres()*TSonataMath::FGauss();
				  }
	       r = gRandom->Rndm();
	       if(r<=par->GetEff()) {
//                First, create a new dighit element
	          dighit_p = new BrDigHit();
	          DigitizedHits->Add(dighit_p);

	          dighit_p->SetID(DigitizedHits->GetEntries());
	          dighit_p->SetIdet(idet);	//+1??to be like FORTRAN
	          dighit_p->SetImod(iplane); //should it be +1 to be like FORTRAN
	          dighit_p->SetWire((Int_t)uwire);
	          dighit_p->SetTime((Int_t)(dU / par->GetDriftv()));
	          dighit_p->SetWidth((Int_t)(par->GetTwoPar()/par->GetDriftv()));
	          dighit_p->SetIfuse(0);
	          dighit_p->SetUsed(0);
	          dighit_p->SetDcpl(iplane);
	          }
	       }
	    }
         }
      }
   }
   }//idet
SortDighit(DigitizedHits);

for(i=0;i<3;i++) fNumGeantHits[i] = GeantHits[i]->GetEntries();
fNumDigHits = DigitizedHits->GetEntries();
if(DebugLevel()>4) {
   for(i=0;i<fNumDigHits;i++) {
      dighit_p = (BrDigHit*)DigitizedHits->At(i);
	  cout<<dighit_p;
      }
   }
}

 void BrDigitizeDC::Decode(BrEventNode* Table)
{

static const float csgn[2] = {(float)-1.0,(float)1.0};
int i,ihit;
float xwire,xdrift;

Int_t NumDigHits,idet;

BrDCPlane* dcpl;
BrDetectorParamsDC *par;
BrDetectorVolume *vol;

BrDataTable* DigitizedHits;
BrDigHit *dighit_p;
BrDCHitPosition *hitpos_p;

//Get stuff from the input table
Char_t TableName[32];
sprintf(TableName,"%s DigitizedDCx0",GetName());
DigitizedHits = Table->GetDataTable(TableName);
NumDigHits = DigitizedHits->GetEntries();
// GetObjectList returns TObjectArray

//Prepare the output table
sprintf(TableName,"%s HitPosDCx0",GetName());
BrDataTable *HitPosition = new BrDataTable(TableName);
Table->AddObject(HitPosition);

par = GetDetectorParamsDC();
if(par == NULL){
   cout << GetName() << "::Digitize; Parameters must be set " << endl;
   return ;
   }

vol = GetDetectorVolume(GetName());

for(idet=0;idet<3;idet++) {
for(ihit=0;ihit<NumDigHits;ihit++) {
   dighit_p = (BrDigHit*)DigitizedHits->At(ihit);
   if( dighit_p->GetIdet() == idet ) {
      if( dighit_p->GetIfuse() == 0 ) {
         Int_t iplane = dighit_p->GetDcpl();
		 dcpl   = par->GetPlaneAt(iplane);
         if( dcpl->GetNwire() == 1 ) xwire = dcpl->GetWirdis() *  (float)dighit_p->GetWire();
         else                        xwire = dcpl->GetWirdis() * ((float)dighit_p->GetWire() + (float)0.5);
         xdrift = par->GetDriftv() * (float)dighit_p->GetTime();
         for( i=0;i<2;i++ ) {
	    hitpos_p = new BrDCHitPosition();
	    HitPosition->Add(hitpos_p);
	    hitpos_p->SetXwire(xwire);
	    hitpos_p->SetXdrift(xdrift);
	    hitpos_p->SetSign(csgn[i]);
	    hitpos_p->SetUhit(xwire + csgn[i]*xdrift);
	    dighit_p->AddHitPos(hitpos_p);
	    if(DebugLevel()>4) cout<<hitpos_p;
	    }
         }
      }
   }
   }//idet
}

 void BrDigitizeDC::SortDighit(BrDataTable* DigitizedHits)
{
Int_t idet,jhit,khit;
Int_t NumDigHits;

BrDigHit* dighit_p1;
BrDigHit* dighit_p2;

NumDigHits = DigitizedHits->GetEntries();

DigitizedHits->Sort(NumDigHits);

idet = 0;		//Place to stop debugger
for( jhit = 0;jhit < NumDigHits-1;jhit++) {
   dighit_p1 = (BrDigHit*)DigitizedHits->At(jhit);
   if( dighit_p1->GetIfuse() == 0 ) {
      for( khit = jhit+1;khit < NumDigHits;khit++ ) {
         dighit_p2 = (BrDigHit*)DigitizedHits->At(khit);
	 if( dighit_p2->GetIfuse() == 0 ) {
	    if( dighit_p1->GetIdet() == dighit_p2->GetIdet() &&
		dighit_p1->GetImod() == dighit_p2->GetImod() &&
		dighit_p1->GetWire() == dighit_p2->GetWire() ) {
	       if( dighit_p1->GetTime()+dighit_p1->GetWidth() > dighit_p2->GetTime() ) {
	          dighit_p1->SetWidth(dighit_p2->GetTime() + dighit_p2->GetWidth() - dighit_p1->GetTime());
		  dighit_p2->SetIfuse(dighit_p1->GetID());
		  }
	       }
	    else if( dighit_p1->GetIdet() == dighit_p2->GetIdet() &&
	       dighit_p1->GetImod() == dighit_p2->GetImod() &&
	       dighit_p1->GetWire() <  dighit_p2->GetWire() ) {
	       //Flemmings GOTO 11 how do we do that here?
	       break;	//This seems to work
	       }
	    else if( dighit_p1->GetIdet() == dighit_p2->GetIdet() &&
		     dighit_p1->GetImod() <  dighit_p2->GetImod() ) {
	       //Flemmings GOTO 11 how do we do that here?
	       break;	//This seems to work
	       }
	    }	        //dighit_p2->ifuse==0?
	 }		//loop over khit
      }			//dighit_p1->ifuse==0?
   }			//loop over jhit
}