//* $Id: */

// -------------------------------------------------------------------------
// -----                    CbmStsDigitize source file             -----
// -----                  Created 08/07/2008  by R. Karabowicz         -----
// -------------------------------------------------------------------------




// Includes from ROOT
#include "TClonesArray.h"
#include "TObjArray.h"
#include "TMath.h"
#include "TF1.h"
#include "TRandom3.h"

#include "TGeoManager.h"
#include "TGeoNode.h"

// Includes from base
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

// Includes from STS
#include "CbmGeoStsPar.h"
#include "CbmStsDigi.h"
#include "CbmStsDigiMatch.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsDigitize.h"
#include "CbmStsPoint.h"
#include "CbmStsSensor.h"
#include "CbmStsSector.h"
#include "CbmStsStation.h"

#include <iostream>
#include <iomanip>
#include <map>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
using std::flush;
using std::pair;
using std::setw;
using std::left;
using std::right;
using std::fixed;
using std::setprecision;
using std::set;
using std::map;
using std::ios_base;
using std::vector;

// -----   Default constructor   ------------------------------------------
CbmStsDigitize::CbmStsDigitize() 
  : FairTask("STS Digitizer", 1),
  fGeoPar(NULL),
  fDigiPar(NULL),
  fPoints(NULL),
  fDigis(NULL),
  fDigiMatches(NULL),
  fRealistic(kFALSE),
  fDigiScheme(NULL),
  fNDigis(0),
  fNMulti(0),
  fNEvents(0),
  fNPoints(0),
  fNOutside(0),
  fNDigisFront(0),
  fNDigisBack(0),
  fStep(0.001),
  fEnergyLossToSignal(0.),
  fFThreshold(8.0),
  fBThreshold(8.0),
  fFNoiseWidth(0.1),
  fBNoiseWidth(0.1),
  fStripDeadTime(10),
  fFNofBits(20),
  fBNofBits(20),
  fFMinStep(0.01),
  fBMinStep(0.01),
  fFNofSteps(0),
  fBNofSteps(0),
  fStripSignalF(NULL),
  fStripSignalB(NULL),
  fTime(0.),
  fTimer(),
  fFChannelPointsMap(),
  fBChannelPointsMap(),
  fPointMap()
{
  fDigiScheme  = new CbmStsDigiScheme();
  Reset();
}
// -------------------------------------------------------------------------

// -----   Standard constructor   ------------------------------------------
CbmStsDigitize::CbmStsDigitize(Int_t iVerbose) 
  : FairTask("STSDigitize", iVerbose), 
  fGeoPar(NULL),
  fDigiPar(NULL),
  fPoints(NULL),
  fDigis(NULL),
  fDigiMatches(NULL),
  fRealistic(kFALSE),
  fDigiScheme(NULL),
  fNDigis(0),
  fNMulti(0),
  fNEvents(0),
  fNPoints(0),
  fNOutside(0),
  fNDigisFront(0),
  fNDigisBack(0),
  fStep(0.001),
  fEnergyLossToSignal(0.),
  fFThreshold(8.0),
  fBThreshold(8.0),
  fFNoiseWidth(0.1),
  fBNoiseWidth(0.1),
  fStripDeadTime(10),
  fFNofBits(20),
  fBNofBits(20),
  fFMinStep(0.01),
  fBMinStep(0.01),
  fFNofSteps(0),
  fBNofSteps(0),
  fStripSignalF(NULL),
  fStripSignalB(NULL),
  fTime(0.),
  fTimer(),
  fFChannelPointsMap(),
  fBChannelPointsMap(),
  fPointMap()
{ 
  fDigiScheme  = new CbmStsDigiScheme();
  Reset();
}
// -------------------------------------------------------------------------

// -----   Constructor with name   -----------------------------------------
CbmStsDigitize::CbmStsDigitize(const char* name, Int_t iVerbose) 
  : FairTask(name, iVerbose), 
  fGeoPar(NULL),
  fDigiPar(NULL),
  fPoints(NULL),
  fDigis(NULL),
  fDigiMatches(NULL),
  fRealistic(kFALSE),
  fDigiScheme(NULL),
  fNDigis(0),
  fNMulti(0),
  fNEvents(0),
  fNPoints(0),
  fNOutside(0),
  fNDigisFront(0),
  fNDigisBack(0),
  fStep(0.001),
  fEnergyLossToSignal(0.),
  fFThreshold(8.0),
  fBThreshold(8.0),
  fFNoiseWidth(0.1),
  fBNoiseWidth(0.1),
  fStripDeadTime(10),
  fFNofBits(20),
  fBNofBits(20),
  fFMinStep(0.01),
  fBMinStep(0.01),
  fFNofSteps(0),
  fBNofSteps(0),
  fStripSignalF(NULL),
  fStripSignalB(NULL),
  fTime(0.),
  fTimer(),
  fFChannelPointsMap(),
  fBChannelPointsMap(),
  fPointMap()
{ 
  fGeoPar      = NULL;
  fDigiPar     = NULL;
  fPoints      = NULL;
  fDigis       = NULL;
  fDigiMatches = NULL;
  fRealistic   = kFALSE;
  fDigiScheme  = new CbmStsDigiScheme();
  Reset();

  fStep = 0.001;

  fFThreshold  = 8.0;
  fBThreshold  = 8.0;
  fFNoiseWidth = 0.1;
  fBNoiseWidth = 0.1;

  fFNofBits    = 20;
  fBNofBits    = 20;
  fFMinStep    = 0.01;
  fBMinStep    = 0.01;
  fStripDeadTime = 10;
  fNEvents = 0.;
}
// -------------------------------------------------------------------------

// -----   Destructor   ----------------------------------------------------
CbmStsDigitize::~CbmStsDigitize() { 
  if ( fGeoPar)    delete fGeoPar;
  if ( fDigiPar)   delete fDigiPar;
  if ( fDigis ) {
    fDigis->Delete();
    delete fDigis;
  }
  if ( fDigiMatches ) {
    fDigiMatches->Delete();
    delete fDigiMatches;
  }
  if ( fDigiScheme ) delete fDigiScheme;
  Reset();
}
// -------------------------------------------------------------------------

// -----   Public method Exec   --------------------------------------------
void CbmStsDigitize::Exec(Option_t* opt) {

//   cout << "threshold = " << fFThreshold << endl;

  // Reset all eventwise counters
  fTimer.Start();
  Reset();
  TF1* digiGausDist = new TF1("digiGausDist","gaus",-5.,5.);

  // Verbose screen output
  if ( fVerbose > 2 ) {
    cout << endl << "-I- " << fName << ": executing event" << endl; 
    cout << "----------------------------------------------" << endl;
  }

  Double_t nPoints  = 0.;
  Double_t nOutside = 0.;
  Double_t nDigisF  = 0.;
  Double_t nDigisB  = 0.;

  // Loop over all StsPoints
  if ( ! fPoints ) {
    cerr << "-W- " << fName << "::Exec: No input array (STSPoint) "
	 << endl;
    cout << "- " << fName << endl;
    return;
  }

  map<CbmStsSensor*, set<Int_t> >::iterator mapIt;
  for (mapIt=fPointMap.begin(); mapIt!=fPointMap.end(); mapIt++)
    ((*mapIt).second).clear();
  
  for (Int_t iPoint=0; iPoint<fPoints->GetEntriesFast(); iPoint++) {
    CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPoint);

    Double_t xin = point->GetXIn();
    Double_t yin = point->GetYIn();
    Double_t zin = point->GetZIn();
    gGeoManager->FindNode(xin,yin,zin);
    TGeoNode* curNode = gGeoManager->GetCurrentNode();
    
    CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curNode->GetName());

    if ( fPointMap.find(sensor) == fPointMap.end() ) {
      cerr << "-E- " << fName << "::Exec:: sensor " << curNode->GetName()
	   << " not found in digi scheme!" << endl;
      continue;
    }
    fPointMap[sensor].insert(iPoint);     
  }
  
  for (Int_t iStation=fDigiScheme->GetNStations(); iStation > 0 ; ) {
    CbmStsStation* station = fDigiScheme->GetStation(--iStation);
    for (Int_t iSector=station->GetNSectors(); iSector > 0 ; ) {
      CbmStsSector* sector = station->GetSector(--iSector);

      map<Int_t, set<Int_t> >::iterator mapCh;

      for (mapCh=fFChannelPointsMap.begin(); mapCh!=fFChannelPointsMap.end(); mapCh++)
	((*mapCh).second).clear();
      for (mapCh=fBChannelPointsMap.begin(); mapCh!=fBChannelPointsMap.end(); mapCh++)
	((*mapCh).second).clear();

      // simulating detector+cables+electronics noise 
      // should be more sophisticated...
      // the question is: sectorwise or sensorwise???
      Int_t nChannels = sector->GetNChannelsFront();
      for (Int_t iChannel=nChannels ; iChannel > 0 ; ) {
// 	fStripSignalF[--iChannel] = fGen->Landau(.1,.02);
// 	fStripSignalB[  iChannel] = fGen->Landau(.1,.02);
	fStripSignalF[--iChannel] = 0.;
	fStripSignalB[  iChannel] = 0.;
// 	fStripSignalF[--iChannel] = TMath::Abs(gRandom->Gaus(0.,fFNoiseWidth));
// 	fStripSignalB[  iChannel] = TMath::Abs(gRandom->Gaus(0.,fBNoiseWidth));
      }
      
      for (Int_t iSensor=sector->GetNSensors(); iSensor > 0 ; ) {
	CbmStsSensor* sensor = sector->GetSensor(--iSensor);
	
	ProduceHitResponse(sensor);
      }

      Int_t   stationNr = sector->GetStationNr();
      Int_t    sectorNr = sector->GetSectorNr();
      Int_t sectorDetId = sector->GetDetectorId();
      
      for ( Int_t ifstr = 0 ; ifstr < nChannels ; ifstr++ ) {
	if ( fStripSignalF[ifstr] < fFThreshold ) continue;
        
	Double_t generator;
	generator = gRandom->Rndm()*100.;	
// 	cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr 
// 	     << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;        
// 	if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ifstr/125])
// 	{
// 	cout << "OCCUPANCYF [" << iStation+1 << "][" << iSector+1 << "][" << ifstr/125 << "] "<< fStripDeadTime*occupancy [iStation][iSector][ifstr/125] << "%" << " generator = "<<generator<< endl;
// 	}
// 	if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ifstr/125]) ifstr = ifstr+2;
	if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ifstr/125]) continue;
	
	Int_t digiFSignal = 1+(Int_t)((fStripSignalF[ifstr]-fFThreshold)/fFMinStep);
	if ( digiFSignal >= fFNofSteps ) digiFSignal = fFNofSteps-1;

// 	new ((      *fDigis)[fNDigis]) CbmStsDigi(stationNr, sectorNr,
// 						  0, ifstr, 
// 						  digiFSignal, 0);
	//						  fStripSignalF[ifstr], 0);
	set<Int_t>::iterator it1;
	set<Int_t> chPnt = fFChannelPointsMap[ifstr];
	Int_t pnt;
	CbmStsDigiMatch* match;
        CbmStsPoint* point = NULL;
        CbmStsDigi* digi;
	if ( chPnt.size() == 0 ) {
// 	  cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr 
// 	       << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
	  new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
	}
	else {
	  for (it1=chPnt.begin(); it1!=chPnt.end(); it1++) {
	    pnt = (*it1);
	    point  = (CbmStsPoint*) fPoints->At(pnt);

 	    Double_t xin = point->GetXIn();
  	    Double_t yin = point->GetYIn();
 	    Double_t zin = point->GetZIn();
    
 	    Double_t xvec = point->GetXOut()-xin;
       	    Double_t yvec = point->GetYOut()-yin;
 	    Double_t zvec = point->GetZOut()-zin;
            Int_t MynofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
    
            Int_t MystepEL = (Int_t) (fEnergyLossToSignal*point->GetEnergyLoss()/(MynofSteps+1));
            xvec = xvec/((Double_t)MynofSteps);  
            yvec = yvec/((Double_t)MynofSteps);  
            zvec = zvec/((Double_t)MynofSteps);  

            for ( Int_t istep = 0 ; istep <= MynofSteps ; istep++ ) {
    	      xin+=xvec;
              yin+=yvec;
       	      zin+=zvec;
   	    }

	    if ( it1==chPnt.begin() ) {
	      match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
	      digi  = new ((      *fDigis)[fNDigis]) CbmStsDigi(pnt, stationNr, sectorNr,
						  0, ifstr, 
						  digiFSignal, 0);
	    }
	    else {
	      match->AddPoint(pnt);
	      fNMulti++;
	      digi->AddIndex(pnt,MystepEL);
	    }
	  }
	}
	fNDigis++;
      }
      
      for ( Int_t ibstr = 0 ; ibstr < nChannels ; ibstr++ ) {
	if ( fStripSignalB[ibstr] < fBThreshold ) continue;

	Double_t generator;
	generator = gRandom->Rndm()*100.;
// 	if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ibstr/125])
// 	{
// 	cout << "OCCUPANCYB [" << iStation+1 << "][" << iSector+1 << "][" << ibstr/125 << "] "<< fStripDeadTime*occupancy [iStation][iSector][ibstr/125] << "%  generator = "<<generator<< endl;
// 	}
// 	if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ibstr/125]) ibstr = ibstr+2;
	if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ibstr/125]) continue;

	Int_t digiBSignal = 1+(Int_t)((fStripSignalB[ibstr]-fBThreshold)/fBMinStep);
	if ( digiBSignal >= fBNofSteps ) digiBSignal = fBNofSteps-1;
// 	new ((      *fDigis)[fNDigis]) CbmStsDigi(stationNr, sectorNr,
// 						  1, ibstr, 
// 						  digiBSignal, 0);
	//						  fStripSignalB[ibstr], 0);
	set<Int_t>::iterator it1;
	set<Int_t> chPnt = fBChannelPointsMap[ibstr];
	Int_t pnt;
        CbmStsPoint* point = NULL;
	CbmStsDigiMatch* match;
        CbmStsDigi* digi;
	if ( chPnt.size() == 0 ) {
// 	  cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr 
// 	       << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
	  new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
	}
	else {
	  for (it1=chPnt.begin(); it1!=chPnt.end(); it1++) {
	    pnt = (*it1);
            point  = (CbmStsPoint*) fPoints->At(pnt);

 	    Double_t xin = point->GetXIn();
  	    Double_t yin = point->GetYIn();
 	    Double_t zin = point->GetZIn();
    
 	    Double_t xvec = point->GetXOut()-xin;
       	    Double_t yvec = point->GetYOut()-yin;
 	    Double_t zvec = point->GetZOut()-zin;
            Int_t MynofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
    
            Int_t MystepEL = (Int_t) (fEnergyLossToSignal*point->GetEnergyLoss()/(MynofSteps+1));
            xvec = xvec/((Double_t)MynofSteps);  
            yvec = yvec/((Double_t)MynofSteps);  
            zvec = zvec/((Double_t)MynofSteps);  

            for ( Int_t istep = 0 ; istep <= MynofSteps ; istep++ ) {
    	      xin+=xvec;
              yin+=yvec;
       	      zin+=zvec;
   	    }


	    if ( it1==chPnt.begin() ) {
	      match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
	      digi  = new ((      *fDigis)[fNDigis]) CbmStsDigi(pnt, stationNr, sectorNr,
						  1, ibstr, 
						  digiBSignal, 0);
	    }
	    else {
	      match->AddPoint(pnt);
	      fNMulti++;
	      digi->AddIndex(pnt,MystepEL);
	    }
	  }
	}
	fNDigis++;
      }
    }
  }
  
  fTimer.Stop();
  cout << "+ " << flush;
  cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8)
       << fixed << right << fTimer.RealTime() 
       << " s, digis " << nDigisF << " / " << nDigisB << endl;
  
  fNEvents     += 1.;
  fNPoints     += Double_t(nPoints);
  fNOutside    += Double_t(nOutside);
  fNDigisFront += Double_t(nDigisF);
  fNDigisBack  += Double_t(nDigisB);
  fTime        += fTimer.RealTime();

}
// -------------------------------------------------------------------------

// -----   Private method ProduceHitResponse   --------------------------------
void CbmStsDigitize::ProduceHitResponse(CbmStsSensor* sensor) {
  set <Int_t> pSet;
  if ( fPointMap.find(sensor) == fPointMap.end() ) {
    cerr << "-E- " << fName << "::Exec:: sensor" 
	 << " not found in digi scheme!" << endl;
    return;
  }
  pSet = fPointMap[sensor];

  Int_t       iPoint = -1;
  CbmStsPoint* point = NULL;

  set<Int_t>::iterator it1;

  for (it1=pSet.begin(); it1!=pSet.end(); it1++) {
    iPoint = (*it1);
    point  = (CbmStsPoint*) fPoints->At(iPoint);

    Double_t xin = point->GetXIn();
    Double_t yin = point->GetYIn();
    Double_t zin = point->GetZIn();
    
    Double_t xvec = point->GetXOut()-xin;
    Double_t yvec = point->GetYOut()-yin;
    Double_t zvec = point->GetZOut()-zin;

    Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
    
    Double_t stepEL = fEnergyLossToSignal*point->GetEnergyLoss()/(nofSteps+1);
    
    xvec = xvec/((Double_t)nofSteps);  
    yvec = yvec/((Double_t)nofSteps);  
    zvec = zvec/((Double_t)nofSteps);  

    for ( Int_t istep = 0 ; istep <= nofSteps ; istep++ ) {
      //      Float_t iFChan = sensor->GetChannelPlus(xin, yin, 0);
      Int_t   iIChan = sensor->GetFrontChannel(xin,yin,zin);
	//(Int_t)iFChan;
      
      if ( iIChan != -1 ) {
	fStripSignalF[iIChan] += stepEL;
	fFChannelPointsMap[iIChan].insert(iPoint);
      }

      //      iFChan = sensor->GetChannelPlus(xin, yin, 1);
      //      iIChan = (Int_t)iFChan;
      iIChan = sensor->GetBackChannel (xin,yin,zin);
      
      if ( iIChan != -1 ) {
	fStripSignalB[iIChan] += stepEL;
	fBChannelPointsMap[iIChan].insert(iPoint);
      }
    
      xin+=xvec;
      yin+=yvec;
      zin+=zvec;
    }
    
  }

}
// -------------------------------------------------------------------------

// -----   Private method FindFiredStrips   --------------------------------
void CbmStsDigitize::FindFiredStrips(CbmStsPoint* pnt,Int_t& nofStr,Int_t*& strips,Double_t*& signals, Int_t side) {

  nofStr = 0;

  Double_t xin = pnt->GetXIn();
  Double_t yin = pnt->GetYIn();
  Double_t zin = pnt->GetZIn();

  gGeoManager->FindNode(xin,yin,zin);
  TGeoNode* curNode = gGeoManager->GetCurrentNode();
  
  CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curNode->GetName());

  Double_t xvec = pnt->GetXOut()-xin;
  Double_t yvec = pnt->GetYOut()-yin;
  Double_t zvec = pnt->GetZOut()-zin;

  Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);

  Double_t stepEL = fEnergyLossToSignal*pnt->GetEnergyLoss()/(nofSteps+1);

  xvec = xvec/((Double_t)nofSteps);  
  yvec = yvec/((Double_t)nofSteps);  
  zvec = zvec/((Double_t)nofSteps);  

  for ( Int_t istep = 0 ; istep <= nofSteps ; istep++ ) {
    Float_t iFChan = sensor->GetChannelPlus(xin, yin, side);
    Int_t   iIChan = (Int_t)iFChan;

    xin+=xvec;
    yin+=yvec;
    zin+=zvec;

    if ( iIChan == -1 ) continue;

    for ( Int_t ifstr = 0 ; ifstr < nofStr ; ifstr++ ) {
      if ( strips[ifstr] == iIChan ) {
	signals[ifstr] += stepEL;
	iIChan = -1;
	break;
      }
    }
    if ( iIChan == -1 ) continue;

    strips [nofStr] = iIChan;
    signals[nofStr] = stepEL;

    nofStr++;

  }
}
// -------------------------------------------------------------------------
 
// -----   Private method SetParContainers   -------------------------------
void CbmStsDigitize::SetParContainers() {

  // Get run and runtime database
  FairRunAna* run = FairRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

  // Get STS geometry parameter container
  fGeoPar = (CbmGeoStsPar*) db->getContainer("CbmGeoStsPar");

  // Get STS digitisation parameter container
  fDigiPar = (CbmStsDigiPar*) db->getContainer("CbmStsDigiPar");
 
}
// -------------------------------------------------------------------------



// -----   Private method Init   -------------------------------------------
InitStatus CbmStsDigitize::Init() {

  // Get input array 
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ) Fatal("Init", "No FairRootManager");
  fPoints = (TClonesArray*) ioman->GetObject("StsPoint");

  // Register output array StsDigi
  fDigis = new TClonesArray("CbmStsDigi",1000);
  ioman->Register("StsDigi", "Digital response in STS", fDigis, kTRUE);

  // Register output array StsDigiMatches
  fDigiMatches = new TClonesArray("CbmStsDigiMatch",1000);
  ioman->Register("StsDigiMatch", "Digi Match in STS", fDigiMatches, kTRUE);

//   fGen = new TRandom3();
//   time_t curtime;
//   time(&curtime);
//   fGen->SetSeed(curtime);

  fStripSignalF = new Double_t[2000];
  fStripSignalB = new Double_t[2000];

  fEnergyLossToSignal    = 200000.;

  fFNofSteps = (Int_t)TMath::Power(2,(Double_t)fFNofBits);
  fBNofSteps = (Int_t)TMath::Power(2,(Double_t)fBNofBits);

  // Build digitisation scheme
  if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
    MakeSets();

    if      (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
    else if (fVerbose >  2) fDigiScheme->Print(kTRUE);
    cout << "-I- " << fName << "::Init: "
	 << "STS digitisation scheme succesfully initialised" << endl;
    cout << "    Stations: " << fDigiScheme->GetNStations() 
	 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: " 
	 << fDigiScheme->GetNChannels() << endl;
    return kSUCCESS;
  }

  return kERROR;

}
// -------------------------------------------------------------------------



// -----   Private method ReInit   -----------------------------------------
InitStatus CbmStsDigitize::ReInit() {

  // Clear digitisation scheme
  fDigiScheme->Clear();

  // Build new digitisation scheme
  if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
    MakeSets();
    return kSUCCESS;
  }

  return kERROR;

}
// -------------------------------------------------------------------------


// -----   Private method MakeSets   ---------------------------------------
void CbmStsDigitize::MakeSets1() {

  fPointMap.clear();
  Int_t nStations = fDigiScheme->GetNStations();
  Double_t fSectorWidth = 0.;
  for (Int_t iStation=0; iStation<nStations; iStation++) {
    CbmStsStation* station = fDigiScheme->GetStation(iStation);
    Int_t nSectors = station->GetNSectors();
    for (Int_t iSector=0; iSector<nSectors; iSector++) {
      CbmStsSector* sector = station->GetSector(iSector);
      Int_t nSensors = sector->GetNSensors();
      for (Int_t iSensor=0; iSensor<nSensors; iSensor++) {
	CbmStsSensor* sensor = sector->GetSensor(iSensor);
	set<Int_t> a;
	fPointMap[sensor] = a;
        fSectorWidth = 10.*sensor->GetLx();
      
        Int_t nofChips = (Int_t)(TMath::Ceil(fSectorWidth/7.5));  // fwidth in mm, 7.5mm = 125(channels)*60mum(pitch)
        Int_t lastChip = (Int_t)(TMath::Ceil(10.*fSectorWidth));
        lastChip = lastChip%75;
        lastChip = (Int_t)(lastChip/.6);
         //     cout << nofChips << " chips on " << iStation+1 << " " << iSector+1 << endl;
        TString addInfo = "";
	if ( nofChips != 8 ) {
	addInfo = Form(", only %d strips",lastChip);
	//	cout << fSectorWidth << " -> " << addInfo.Data() << endl;
        }
   
      
        for ( Int_t iChip = 0 ; iChip < nofChips ; iChip++ ) {
          occupancy [iStation][iSector][iChip] = 3.;
//           cout << "OCCUPANCY [" << iStation+1 << "][" << iSector+1 << "][" << iChip << "] "<< occupancy [iStation][iSector][iChip] << "%" << endl;
        }
      }
    }
  }
  fFChannelPointsMap.clear();
  fBChannelPointsMap.clear();
  for ( Int_t ichan = 2000 ; ichan > 0 ; ) {
    set<Int_t> a;
    fFChannelPointsMap[--ichan] = a;
    set<Int_t> b;
    fBChannelPointsMap[  ichan] = b;
  }
}
// -------------------------------------------------------------------------
void CbmStsDigitize::MakeSets() {


  fPointMap.clear();
  Int_t nStations = fDigiScheme->GetNStations();
      
  TH1F* fhFNofDigisPChip[10][1000][20];
  TH1F* fhBNofDigisPChip[10][1000][20];
  TString qaFileName;
  qaFileName = "occup.sts.reco.root";
  cout << "Occupancy read from file: \"" << qaFileName.Data() << "\"" << endl;
  TFile* occuF = TFile::Open(qaFileName.Data());  
  if ( !occuF ) {
    cout << "sorry, no file" << endl; 
    MakeSets1();
    return; 
  }
  TString directoryName = "STSFindHitsQA";
  
  Double_t fSectorWidth = 0.;
  
  for (Int_t iStation=0; iStation<nStations; iStation++) {
    CbmStsStation* station = fDigiScheme->GetStation(iStation);
    Int_t nSectors = station->GetNSectors();
    
    for (Int_t iSector=0; iSector<nSectors; iSector++) {
      CbmStsSector* sector = station->GetSector(iSector);
      Int_t nSensors = sector->GetNSensors();
      Int_t nChannels = sector->GetNChannelsFront();
      
      for (Int_t iSensor=0; iSensor<nSensors; iSensor++) {
	CbmStsSensor* sensor = sector->GetSensor(iSensor);
	set<Int_t> a;
	fPointMap[sensor] = a;
	fSectorWidth = 10.*sensor->GetLx();
      
      
        Int_t nofChips = (Int_t)(TMath::Ceil(fSectorWidth/7.5));  // fwidth in mm, 7.5mm = 125(channels)*60mum(pitch)
        Int_t lastChip = (Int_t)(TMath::Ceil(10.*fSectorWidth));
        lastChip = lastChip%75;
        lastChip = (Int_t)(lastChip/.6);
         //     cout << nofChips << " chips on " << iStation+1 << " " << iSector+1 << endl;
        TString addInfo = "";
	if ( nofChips != 8 ) {
	addInfo = Form(", only %d strips",lastChip);
	//	cout << fSectorWidth << " -> " << addInfo.Data() << endl;
        }
   
      
        for ( Int_t iChip = 0 ; iChip < nofChips ; iChip++ ) {
          fhFNofDigisPChip[iStation][iSector][iChip]=(TH1F*)occuF->Get(Form("%s/Station%d/hNofFiredDigisFSt%dSect%dChip%d",directoryName.Data(),iStation+1,iStation+1,iSector+1,iChip+1));
          fhBNofDigisPChip[iStation][iSector][iChip]=(TH1F*)occuF->Get(Form("%s/Station%d/hNofFiredDigisBSt%dSect%dChip%d",directoryName.Data(),iStation+1,iStation+1,iSector+1,iChip+1));
          occupancy [iStation][iSector][iChip] = 100.*fhFNofDigisPChip[iStation][iSector][iChip]->GetMean()/125.;
          occupancy [iStation][iSector][iChip] = 100.*fhBNofDigisPChip[iStation][iSector][iChip]->GetMean()/125.;
// 	    if ( !occuF ) {
// 	      occupancy [iStation][iSector][iChip] = 3.;
// 	    }
//           cout << "OCCUPANCY [" << iStation+1 << "][" << iSector+1 << "][" << iChip << "] "<< occupancy [iStation][iSector][iChip] << "%" << endl;
        }
        

      }
    }
  }
  fFChannelPointsMap.clear();
  fBChannelPointsMap.clear();
  for ( Int_t ichan = 2000 ; ichan > 0 ; ) {
    set<Int_t> a;
    fFChannelPointsMap[--ichan] = a;
    set<Int_t> b;
    fBChannelPointsMap[  ichan] = b;
  }
 
}
// -------------------------------------------------------------------------

// -----   Private method Reset   ------------------------------------------
void CbmStsDigitize::Reset() {
  //  fNPoints = fNOutside = fNDigisFront = fNDigisBack = fTime = 0.;
  fNDigis = fNMulti = 0;
  fFChannelPointsMap.clear();
  fBChannelPointsMap.clear();
  //  if ( fDigis ) fDigis->Clear();
  //  if ( fDigiMatches ) fDigiMatches->Clear();
  if ( fDigis ) fDigis->Delete();
  if ( fDigiMatches ) fDigiMatches->Delete();
}
// -------------------------------------------------------------------------


// -----   Virtual method Finish   -----------------------------------------
void CbmStsDigitize::Finish() {
  cout << endl;
  cout << "============================================================"
       << endl;
  cout << "===== " << fName << ": Run summary " << endl;
  cout << "===== " << endl;
  cout << "===== Events processed          : " << setw(8) << fNEvents << endl;
  cout.setf(ios_base::fixed, ios_base::floatfield);
  cout << "===== Real time per event       : " 
       << setw(8) << setprecision(4) 
       << fTime / fNEvents << " s" << endl;
  cout << "===== StsPoints per event       : " 
       << setw(8) << setprecision(2) 
       << fNPoints / fNEvents << endl;
  cout << "===== Outside hits per event    : " 
       << setw(8) << setprecision(2) 
       << fNOutside / fNEvents << " = " 
       << setw(6) << setprecision(2) 
       << fNOutside / fNPoints * 100. << " %" << endl;
  cout << "===== Front digis per point     : " 
       << setw(8) << setprecision(2) 
       << fNDigisFront / (fNPoints-fNOutside) << endl;
  cout << "===== Back digis per point      : " 
       << setw(8) << setprecision(2) 
       << fNDigisBack / (fNPoints-fNOutside) << endl;
  cout << "============================================================"
       << endl;
	
}					       
// -------------------------------------------------------------------------



ClassImp(CbmStsDigitize)