// -------------------------------------------------------------------------
// -----                PndMvdNoiseProducer source file                -----
// -----                  Created 01.07.08  by R.Kliemt                -----
// -------------------------------------------------------------------------

//#include <cmath>

#include "TClonesArray.h"
#include "TGeoNode.h"

#include "FairRootManager.h"
#include "FairGeoVolume.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "FairGeoNode.h"

#include "PndMvdNoiseProducer.h"
#include "PndSdsMCPoint.h"
#include "PndSdsDigiStrip.h"
#include "PndSdsDigiPixel.h"
#include "PndStringSeparator.h"

#include "PndSdsIdealChargeConversion.h"
#include "PndSdsTotChargeConversion.h"

#include "PndDetectorList.h"

// -----   Default constructor   -------------------------------------------
PndMvdNoiseProducer::PndMvdNoiseProducer() :
FairTask("Charge Noise Producer"), fGeoH(0)
{
  fBranchName 	= "MVDStripDigis";
  fPersistance = kTRUE;
  //  fGeoH = NULL;

  fTimeOrderedDigi = kFALSE;
  fBranchName = "";
  fDigiStripArray = NULL;
  fDigiPixelArray = NULL;
  fDigiPixelBuffer = NULL;
  fDigiStripBuffer = NULL;
  fDigiParRect = NULL;
  fDigiParTrap = NULL;
  fDigiParPix = NULL;
  fTotDigiParRect = NULL;
  fTotDigiParTrap = NULL;
  fTotDigiParPix = NULL;
  fMCEventheader = NULL;
  fStripRectChargeConv = NULL;
  fStripTrapChargeConv = NULL;
  fCurrentChargeConv = NULL;
  fPixChargeConv = NULL;
  fNoiseSpread = 0;
  fThreshold = 0;
  fPreviousTime = 0.;
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
PndMvdNoiseProducer::~PndMvdNoiseProducer()
{
	if (fStripRectChargeConv) delete fStripRectChargeConv;
	if (fStripTrapChargeConv) delete fStripTrapChargeConv;
	if (fPixChargeConv) delete fPixChargeConv;
}

// -----   Public method Init   --------------------------------------------
InitStatus PndMvdNoiseProducer::Init()
{
 // Get RootManager
  FairRootManager* ioman = FairRootManager::Instance();
  
  if ( ! ioman )
  {
	std::cout << " -E- PndMvdNoiseProducer::Init: RootManager not instantiated!" << std::endl;
	return kFATAL;
  }
  
  // Get input array

  fDigiStripBuffer = (PndSdsDigiStripWriteoutBuffer*)FairRootManager::Instance()->RegisterWriteoutBuffer("MVDStripDigis", new PndSdsDigiStripWriteoutBuffer("MVDStripDigis", "MVD", kTRUE));
  
//  fDigiPixelArray = (TClonesArray*) ioman->GetObject("MVDPixelDigis");
  fDigiPixelBuffer = new PndSdsDigiPixelWriteoutBuffer("MVDPixelDigis", "MVD", kTRUE);
  fDigiPixelBuffer = (PndSdsDigiPixelWriteoutBuffer*)FairRootManager::Instance()->RegisterWriteoutBuffer("MVDPixelDigis", fDigiPixelBuffer);
  fDigiPixelBuffer->ActivateBuffering(fTimeOrderedDigi);

  
  fMCEventheader = (FairMCEventHeader*) ioman->GetObject("MCEventHeader.");
  if ( ! fMCEventheader ){
	Warning("Init","Did not find the MC event header, assume 50ns of noise clockticks per call of Exec().");
  }
  fPreviousTime=0.;
  
  FillSensorLists();
  
//  fFEModel = new PndSdsFESimple();

  if(fVerbose>0)
  {
	std::cout <<" -I- PndMvdNoiseProducer: Registered Sensors: "
	<<fStripRectLIds.size()+fStripRectSIds.size()<<"xStripRect "
	<<fStripTrapIds.size()<<"xStripTrap "
	<<fPixelIds2.size()<<"xPixel"
	<<std::endl;
  }
  std::cout << " -I- PndMvdNoiseProducer: Intialisation successfull" << std::endl;
  
  if (fDigiParRect->GetChargeConvMethod() == 0){
	if(fVerbose>0) Info("Init()","ideal charge conversion for rect. strips");
	fStripRectChargeConv = new PndSdsIdealChargeConversion(fDigiParRect->GetNoise());
  }
  else if (fDigiParRect->GetChargeConvMethod() == 1){
	if(fVerbose>0) Info("Init()","use TOT charge conversion for rect. strips");
	fStripRectChargeConv = new PndSdsTotChargeConversion(
							     fTotDigiParRect->GetChargingTime(),
							     fTotDigiParRect->GetConstCurrent(),
							     fDigiParRect->GetThreshold(),
							     fTotDigiParRect->GetClockFrequency(),
							     fVerbose);
  }
  else Fatal ("Init()","rect. strips: charge conversion method not defined!");
  
  if (fDigiParTrap->GetChargeConvMethod() == 0){
    if(fVerbose>0) Info("Init()","ideal charge conversion for trap. strips");
    fStripTrapChargeConv = new PndSdsIdealChargeConversion(fDigiParTrap->GetNoise());
  }
  else if (fDigiParTrap->GetChargeConvMethod() == 1){
    if(fVerbose>0) Info("Init()","use TOT charge conversion for trap. strips");
    fStripTrapChargeConv = new PndSdsTotChargeConversion(
							 fTotDigiParTrap->GetChargingTime(),
							 fTotDigiParTrap->GetConstCurrent(),
							 fDigiParTrap->GetThreshold(),
							 fTotDigiParTrap->GetClockFrequency(),
							 fVerbose);
  }
  else Fatal ("Init()","trap. strips: charge conversion method not defined!");
  
  if (fDigiParPix->GetChargeConvMethod() == 0){
	if(fVerbose>0) Info("Init()","ideal charge conversion for pixel part");
	fPixChargeConv = new PndSdsIdealChargeConversion(fDigiParPix->GetNoise());
  }
  else if (fDigiParPix->GetChargeConvMethod() == 1){
    if(fVerbose>0) Info("Init()","use TOT charge conversion for pixel part");
    fPixChargeConv = new PndSdsTotChargeConversion(
						   fTotDigiParPix->GetChargingTime(),
						   fTotDigiParPix->GetConstCurrent(),
						   fDigiParPix->GetThreshold(),
						   fTotDigiParPix->GetClockFrequency(),
						   fVerbose);
  }
  else Fatal ("Init()","pixel part: charge conversion method not defined!");
  
  return kSUCCESS;
}

void PndMvdNoiseProducer::FillSensorLists()
{
  TObjArray* sensorNames = fGeoH->GetSensorNames();

  for (int i = 0; i < sensorNames->GetEntries(); i++)
  {
		TString volpath = ((TObjString*)(sensorNames->At(i)))->GetString();
    if(!volpath.Contains("Mvd")) continue;
		PndStringSeparator sep(volpath.Data(), "/");
		std::vector<std::string> volvec = sep.GetStringVector();
		TString volname = volvec[volvec.size()-1].c_str();
		if(fVerbose>2)std::cout << "VolName: " << volname.Data();
		if(volname.Contains("Active"))
		{
		  if(volname.Contains("RectL")) {fStripRectLIds.push_back(i); if(fVerbose>2)std::cout << " \tAdded to StripRectL" << std::endl;}
		  if(volname.Contains("RectS")) {fStripRectSIds.push_back(i); if(fVerbose>2)std::cout << " \tAdded to StripRectS" << std::endl;}
		  if(volname.Contains("Trap")) {fStripTrapIds.push_back(i); if(fVerbose>2)std::cout << " \tAdded to StripTrap" << std::endl;}
		  if(volname.Contains("Pixel"))
		  {
        if(volname.Contains("2")) {fPixelIds2.push_back(i); if(fVerbose>2)std::cout << " \tAdded to Pixel 2" << std::endl;}
        if(volname.Contains("4")) {fPixelIds4.push_back(i); if(fVerbose>2)std::cout << " \tAdded to Pixel 4" << std::endl;}
        if(volname.Contains("5")) {fPixelIds5.push_back(i); if(fVerbose>2)std::cout << " \tAdded to Pixel 5" << std::endl;}
        if(volname.Contains("6")) {fPixelIds6.push_back(i); if(fVerbose>2)std::cout << " \tAdded to Pixel 6" << std::endl;}
		  }
		}
	}
}

// -------------------------------------------------------------------------
void PndMvdNoiseProducer::SetParContainers()
{
  if ( fGeoH == NULL ) {
    fGeoH = PndGeoHandling::Instance();
  }

  fGeoH->SetParContainers();
  // Get Base Container
  FairRun* ana = FairRun::Instance();
  FairRuntimeDb* rtdb=ana->GetRuntimeDb();
  fDigiParRect = (PndSdsStripDigiPar*)(rtdb->getContainer("MVDStripDigiParRect"));
  fDigiParTrap = (PndSdsStripDigiPar*)(rtdb->getContainer("MVDStripDigiParTrap"));
  fDigiParPix  = (PndSdsPixelDigiPar*)(rtdb->getContainer("MVDPixelDigiPar"));
  fTotDigiParRect = (PndSdsTotDigiPar*)(rtdb->getContainer("MVDStripTotDigiParRect"));
  fTotDigiParTrap = (PndSdsTotDigiPar*)(rtdb->getContainer("MVDStripTotDigiParTrap"));
  fTotDigiParPix  = (PndSdsTotDigiPar*)(rtdb->getContainer("MVDPixelTotDigiPar"));
}


// -----   Public method Exec   --------------------------------------------
void PndMvdNoiseProducer::Exec(Option_t* opt)
{
  //   TObjArray* activeSensors = fGeoPar->GetGeoSensitiveNodes();
  Int_t nrCh=0,rnd=0,
  nrFE=0,sens=0,
  nrSensors=0,fe=0,
  chanmax=0,chan=0,
  col=0,row=0,
  chanwhite=0,
  charge=0,
  nNoisyStripRects=0,
  nNoisyStripTraps=0,
  nNoisyPixels=0;
  Double_t xfrac=0.,cycles=1.;
  Int_t did=-1;
  
  // *** Strip Rect Long ***
  // how many channels left?
  nrCh = fDigiParRect->GetNrFECh();
  nrFE = fDigiParRect->GetNrBotFE() + fDigiParRect->GetNrTopFE();
  nrSensors = fStripRectLIds.size();
  chanmax = nrCh * nrFE * nrSensors;
  // Get Number of Channels fired from noise
  xfrac = CalcDistFraction(fDigiParRect->GetNoise(),fDigiParRect->GetThreshold());
  cycles = CalcReadoutCycles(fDigiParRect->GetFeBusClock());
  chanwhite = gRandom->Poisson(xfrac*chanmax);
  if(fVerbose>1) std::cout << " -I- PndMvdNoiseProducer: RECT <N> = " << xfrac*cycles*chanmax
    << " leading to " << chanwhite << " noisy digis of " << chanmax
    << " total channels" << std::endl;
  for(Int_t i = 0;i < chanwhite;i++)
  {
    // randomize the channel numbers & sensors
    rnd = gRandom->Integer(chanmax);
    sens = rnd/(nrFE*nrCh);
    rnd = rnd % (nrFE*nrCh);
    fe = rnd/nrCh; 
    chan = rnd % nrCh;
    // calculate a charge deposit above threshold
    charge = CalcChargeAboveThreshold(fDigiParRect->GetNoise(),fDigiParRect->GetThreshold());
    did = fStripRectLIds[sens];
    fCurrentChargeConv = fStripRectChargeConv;
    AddDigiStrip(nNoisyStripRects,-1,did,fe,chan,charge);
  }
  
  // *** Strip Rect Short ***
  // how many channels left?
  nrCh = fDigiParRect->GetNrFECh();
  nrFE = fDigiParRect->GetNrBotFE() + fDigiParRect->GetNrTopFE()/2 ;
  nrSensors = fStripRectSIds.size();
  chanmax = nrCh * nrFE * nrSensors;
  // Get Number of Channels fired from noise
  xfrac = CalcDistFraction(fDigiParRect->GetNoise(),fDigiParRect->GetThreshold());
  chanwhite = gRandom->Poisson(xfrac*cycles*chanmax);
  if(fVerbose>1) std::cout << " -I- PndMvdNoiseProducer: RECT <N> = " << xfrac*cycles*chanmax
    << " leading to " << chanwhite << " noisy digis of " << chanmax
    << " total channels" << std::endl;
  for(Int_t i = 0;i < chanwhite;i++)
  {
    // randomize the channel numbers & sensors
    rnd = gRandom->Integer(chanmax);
    sens = rnd/(nrFE*nrCh);
    rnd = rnd % (nrFE*nrCh);
    fe = rnd/nrCh; 
    if(fe >= fDigiParRect->GetNrTopFE()/2) {
      fe += fDigiParRect->GetNrTopFE()/2; // shift to avoid the unused fe numbers.
    }
    chan = rnd % nrCh;
    // calculate a charge deposit above threshold
    charge = CalcChargeAboveThreshold(fDigiParRect->GetNoise(),fDigiParRect->GetThreshold());
    did = fStripRectSIds[sens];
    fCurrentChargeConv = fStripRectChargeConv;
    AddDigiStrip(nNoisyStripRects,-1,did,fe,chan,charge);
  }
  
  // *** Strip Trapezoids ***
  nrCh = fDigiParTrap->GetNrFECh();
  nrFE = fDigiParTrap->GetNrBotFE() + fDigiParTrap->GetNrTopFE();
  nrSensors = fStripTrapIds.size();
  chanmax = nrCh * nrFE * nrSensors;
  xfrac = CalcDistFraction(fDigiParTrap->GetNoise(),fDigiParTrap->GetThreshold());
  cycles = CalcReadoutCycles(fDigiParTrap->GetFeBusClock());
  chanwhite = gRandom->Poisson(xfrac*cycles*chanmax);
  if(fVerbose>1) std::cout << " -I- PndMvdNoiseProducer: TRAP <N> = " << xfrac*cycles*chanmax
    << " leading to " << chanwhite << " noisy digis of " << chanmax
    << " total channels" << std::endl;
  for(Int_t i = 0;i < chanwhite;i++)
  {
    rnd = gRandom->Integer(chanmax);
    sens = rnd/(nrFE*nrCh);
    rnd = rnd % (nrFE*nrCh);
    fe = rnd/nrCh;
    chan = rnd % nrCh;
    charge = CalcChargeAboveThreshold(fDigiParTrap->GetNoise(),fDigiParTrap->GetThreshold());
    did = fStripTrapIds[sens];
    fCurrentChargeConv = fStripTrapChargeConv;
    AddDigiStrip(nNoisyStripTraps,-1,did,fe,chan,charge);
  }
  
  // *** Pixel Sensors ***
  nrCh = fDigiParPix->GetFECols()*fDigiParPix->GetFERows();
  Int_t pixx2=2*fPixelIds2.size(), 
  pixx4=4*fPixelIds4.size(), 
  pixx5=5*fPixelIds5.size(), 
  pixx6=6*fPixelIds6.size();
  nrFE = pixx2 + pixx4 + pixx5 + pixx6;
  chanmax = nrCh * nrFE;
  xfrac = CalcDistFraction(fDigiParPix->GetNoise(),fDigiParPix->GetThreshold());
  cycles = CalcReadoutCycles(fDigiParPix->GetFeBusClock());
  chanwhite = gRandom->Poisson(xfrac*cycles*chanmax);
  if(fVerbose>1) std::cout << " -I- PndMvdNoiseProducer: PIXEL <N> = " << xfrac*cycles*chanmax
    << " leading to " << chanwhite << " noisy digis of " << chanmax
    << " total channels" << std::endl;
  for(Int_t i = 0;i < chanwhite;i++)
  {
    charge = CalcChargeAboveThreshold(fDigiParPix->GetNoise(),fDigiParPix->GetThreshold());
    rnd = gRandom->Integer(chanmax);
    chan = rnd%nrCh;
    col = chan%fDigiParPix->GetFECols();
    row = chan/fDigiParPix->GetFECols();
    fe  = rnd/nrCh;
    if(fe >= (Int_t)(pixx2 + pixx4 + pixx5) )
    {
      fe = fe - pixx2 - pixx4 - pixx5;
      sens = fe/6;
      did = fPixelIds6[sens];
      fe = fe%6;
      //if(fe>6) fe=fe-6+10; //0-9 one row of FE, 10-19 2nd row of FE
    } else if( fe >= (Int_t)(pixx2 + pixx4) )
    {
      fe = fe - pixx2 - pixx4;
      sens = fe/5;
      did = fPixelIds5[sens];
      fe = fe%5;
      //if(fe>4) fe=fe-4+10; //0-9 one row of FE, 10-19 2nd row of FE
    } else if( fe >= (Int_t)(pixx2) )
    {
      fe = fe -pixx2;
      sens = fe/4;
      did = fPixelIds4[sens];
      fe = fe%4;
    } else
    {
      sens = fe/2;
      did = fPixelIds2[sens];
      fe = fe%2;
    }
    
    AddDigiPixel(nNoisyPixels,-1,did,fe,col,row,charge);
  }
  
  fPreviousTime = FairRootManager::Instance()->GetEventTime();

  // *** The End ***
  if(fVerbose>0)
  {
    std::cout <<" -I- PndMvdNoiseProducer: Noise produced\t"
    <<nNoisyStripRects <<"xStripRect\t"
    <<nNoisyStripTraps <<"xStripTrap\t"
    <<nNoisyPixels <<"xPixels"<<std::endl;
  }
}
// -------------------------------------------------------------------------

// -------------------------------------------------------------------------
Double_t PndMvdNoiseProducer::CalcDistFraction(Double_t spread, Double_t threshold)
{
  // mean fraction of fireing digis
  return TMath::Erfc( threshold / (TMath::Sqrt2()*spread) );
}

Int_t PndMvdNoiseProducer::CalcChargeAboveThreshold(Double_t spread,Double_t threshold)
{
  // pick any charge from the gaussian tail above threshold
  Double_t temp=0.;
  temp = TMath::Gaus(threshold,0,spread); // maximum gauss at threshold
  temp = gRandom->Uniform(0.,temp);       // random value in y
  temp = -2.*spread*spread*log(temp);     // get x value (recalc charge)
  temp = sqrt(temp);
  return (Int_t)ceil(temp);
}

Double_t PndMvdNoiseProducer::CalcReadoutCycles(Double_t clock)
{ // time [ns], clock [MHz]
  Double_t cycles=1.;
  Double_t timewindow=0.;
  if (clock > 0){
    if (fMCEventheader!=0) {
      timewindow = FairRootManager::Instance()->GetEventTime();
      timewindow -= fPreviousTime;
    } else {
      timewindow = 50.; // how many ns do we suppress readout of zeros?
    }
  }
  if(fVerbose>1) printf(" -I- PndMvdNoiseProducer::CalcReadoutCycles(): %g cycles (%gMHz,%gns)\n",cycles,clock,timewindow);
  return cycles;
}

void PndMvdNoiseProducer::AddDigiStrip(Int_t &noisies, Int_t iPoint, Int_t sensorID, Int_t fe, Int_t chan, Double_t charge)
{
  //Double_t tempcharge = 0.;
  //Bool_t found = kFALSE;
  Int_t detID = -1; // no source mc branch
  //PndSdsDigiStrip* aDigi = 0;
//  FairMCEventHeader* MCevtHeader = (FairMCEventHeader*)FairRootManager::Instance()->GetObject("MCEventHeader.");

//  if (fTimeOrderedDigi == kFALSE){
//	  Int_t iStrip = fDigiStripArray->GetEntriesFast();
//
//	  for(Int_t kstr = 0; kstr < iStrip && found == kFALSE; kstr++)
//	  {
//		aDigi = (PndSdsDigiStrip*)fDigiStripArray->At(kstr);
//		if (aDigi->GetSensorID() == sensorID &&
//			aDigi->GetFE() == fe &&
//			aDigi->GetChannel() == chan )
//		{
//		  tempcharge = fCurrentChargeConv->DigiValueToCharge(*aDigi);
//		  aDigi->SetCharge( fCurrentChargeConv->ChargeToDigiValue(charge + tempcharge) );
//		  aDigi->AddIndex(iPoint);
//		  found = kTRUE;
//		}
//	  }
//	  if(found == kFALSE){
//		  //TODO: get a reasonable timestamp fake for the noise
//		std::vector<Int_t> indices;
//		indices.push_back(iPoint);
//		new ((*fDigiStripArray)[iStrip]) PndSdsDigiStrip(indices,detID,sensorID,fe,chan,fCurrentChargeConv->ChargeToDigiValue(charge), FairRootManager::Instance()->GetEventTime()) ;
//		noisies++;
//		if(fVerbose>2) std::cout
//		  << " -I- PndSdsNoiseProducer: Added StripTrap Digi at: FE=" << fe
//		  << ", channel=" << chan << ", charge=" << charge<< " e"
//		  << ", in sensor \n" << sensorID <<std::endl;
//	  }
//  }
//  else{
		std::vector<Int_t> indices;
		indices.push_back(iPoint);
		PndSdsDigiStrip* tempStrip = new PndSdsDigiStrip(indices,detID,sensorID,fe,chan,fCurrentChargeConv->ChargeToDigiValue(charge), FairRootManager::Instance()->GetEventTime()) ;
		noisies++;
		fDigiStripBuffer->FillNewData(tempStrip, FairRootManager::Instance()->GetEventTime(), FairRootManager::Instance()->GetEventTime() + 10);
		delete(tempStrip);
//  }
}
// -------------------------------------------------------------------------
void PndMvdNoiseProducer::AddDigiPixel(Int_t &noisies, Int_t iPoint, Int_t sensorID, Int_t fe, Int_t col, Int_t row, Double_t charge)
{
  //Double_t tempcharge = 0.;
  //Bool_t found = kFALSE;
  Int_t detID = -1; //no source mc branch

  std::vector<Int_t> indices;
  indices.push_back(iPoint);
  PndSdsDigiPixel* tempPixel = new PndSdsDigiPixel(indices,detID,sensorID,fe,col,row,fPixChargeConv->ChargeToDigiValue(charge), fPixChargeConv->GetTimeStamp(0, charge,FairRootManager::Instance()->GetEventTime()));//FairRootManager::Instance()->GetEventTime()) ;
  if (fPixChargeConv->GetTimeWalk((Int_t)tempPixel->GetCharge()) < 1E5){
		tempPixel->SetTimeStamp(tempPixel->GetTimeStamp() - fPixChargeConv->GetTimeWalk((Int_t)tempPixel->GetCharge()));
		tempPixel->SetTimeStampError(fPixChargeConv->GetTimeStampErrorAfterCorrection());
  }
  fDigiPixelBuffer->FillNewData(tempPixel,fPixChargeConv->ChargeToDigiValue(charge)*6 + FairRootManager::Instance()->GetEventTime(), FairRootManager::Instance()->GetEventTime());

  delete(tempPixel);
  //  std::cout << "DataInBuffer: " << fDigiPixelBuffer->GetNData() << std::endl;
}

void PndMvdNoiseProducer::FinishEvent()
{
  // called after all Tasks did their Exex() and the data is copied to the file

  FinishEvents();
}
// -------------------------------------------------------------------------




ClassImp(PndMvdNoiseProducer)