// -------------------------------------------------------------------------
// -----                    CbmMvdSensorBufferTask  source file                -----
// -------------------------------------------------------------------------


/*
 *
 * ____________________________________________________________________________________________
  
 * ____________________________________________________________________________________________
 * --------------------------------------------------------------------------------------------
 */



// Includes from MVD
#include "CbmMvdSensorBufferTask.h"


#include "CbmMvdGeoPar.h"
#include "CbmMvdHit.h"
#include "CbmMvdHitMatch.h"
#include "CbmMvdPileupManager.h"
#include "CbmMvdPoint.h"
#include "CbmMvdStation.h"
#include "CbmMvdSensor.h"
#include "CbmMvdPointToHitConverter.h"

#include "../buffers/CbmMvdSensorFrameBuffer.h"
#include "../buffers/CbmMvdSensorTrackingBuffer.h"
//#include "omp.h"

// Includes from base
#include "FairGeoNode.h"
#include "FairRootManager.h"
#include "FairRunSim.h"
#include "FairRuntimeDb.h"
#include "CbmMCTrack.h"

// Includes from ROOT
#include "TArrayD.h"
#include "TClonesArray.h"
#include "TGeoManager.h"
//#include "TGeoShape.h"
#include "TGeoTube.h"
#include "TObjArray.h"
#include "TRandom3.h"
#include "TString.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"

// Includes from C++
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>

//#include "gsl/gsl_sf_erf.h"
//#include "CLHEP/Random/RandGauss.h"
//#include "CLHEP/Random/RandPoisson.h"
//#include "CLHEP/Random/RandFlat.h"


using std::cout;
using std::endl;
using std::map;
using std::setw;
using std::left;
using std::right;
using std::fixed;
using std::pair;
using std::setprecision;
using std::ios_base;
using std::vector;


// -----   Default constructor   ------------------------------------------
CbmMvdSensorBufferTask::CbmMvdSensorBufferTask()
    : FairTask("MVDSensorBuffer")
{
  
    bWriteOutput   = false;
    fMode          = 0;
    fBranchName    = "MvdPoint";
    fDigis         = new TClonesArray("CbmMvdDigi");
    fPixelCharge   = new TClonesArray("CbmMvdPixelCharge");
    fPileupManager = NULL;
    fDeltaManager  = NULL;
    fRandGen.SetSeed(2736);
    fEvent       = 0;
    fTime        = 0.;
    fSigmaX      = 0.0005;
    fSigmaY      = 0.0005;
    fReadoutTime = 0.00005;
    fNPileup     = 0;
    fNDeltaElect = 0;
    fBgFileName    = "";
    fDeltaFileName = "";
    fBgBufferSize    = 1000;
    fDeltaBufferSize = 10000;

    fPoints=new TRefArray();
    //fFluctuate = new MyG4UniversalFluctuationForSi();

    fEpiTh         = 0.0014;
    fSegmentLength = 0.0001;
    fDiffusionCoefficient = 0.0055; // correspondes to the sigma of the gauss with the max drift length
    fElectronsPerKeV = 276; //3.62 eV for e-h creation
    fWidthOfCluster  = 3.5; // in sigmas
    fPixelSizeX = 0.0030; // in cm
    fPixelSizeY = 0.0030;
    fCutOnDeltaRays  = 0.00169720;  //MeV
    fChargeThreshold = 150; //electrons change 1 to 150
    fFanoSilicium    = 0.115;
    fEsum            = 0;
    fSegmentDepth    = 0;
    fCurrentTotalCharge      = 0;
    fCurrentParticleMass     = 0;
    fCurrentParticleMomentum = 0;
    fPixelScanAccelerator    = 0;
    fLandauRandom=new TRandom3();

    fShowDebugHistos = kFALSE;

    fPixelSize = 0.0010;
    fPar0 = 3.42157e+02;
    fPar1 = 7.67981e-01;
    fPar2 = 0;

    fLandauMPV=1013.;
    fLandauSigma=159.2;
    fLandauGain=1.56;

    /*fLorentzY0=-21.74022;
    fLorentzXc=0.; // -0.04133
    fLorentzW=1.28456;
    fLorentzA=1008.79877; */




    /*
    //Mimosa17 - line in cluster
    fLorentzY0=-1.7;
    fLorentzXc=0.; // -0.04133
    fLorentzW=1.107;
    fLorentzA=443.35;
    */

    /*
    //Mimosa17 - column in cluster
    fLorentzY0=-0.16;
    fLorentzXc=0.; // -0.04133
    fLorentzW=0.905;
    fLorentzA=357.02;
    */

    /*
    //params adjusted from "intermediate" mpv values for mimosa17
    fLorentzY0=-1.99;
    fLorentzXc=0.;
    fLorentzW=1.056;
    fLorentzA=416.4;
    */

    /*
    //mimosa18 - data from column
    fLorentzY0=-6;
    fLorentzXc=0.;
    fLorentzW=0.95;
    fLorentzA=430.4;
    */

    //mimosa18 - average column + line
    fLorentzY0=-6.1;
    fLorentzXc=0.;
    fLorentzW=1.03;
    fLorentzA=477.2;

    //fLorentzNorm=0.00013010281679422413;
    fLorentzNorm=1;
    
    SetPixelSize(18.4);
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
CbmMvdSensorBufferTask::CbmMvdSensorBufferTask(const char* name, Int_t iMode,
				 Int_t iVerbose)
: FairTask(name, iVerbose)
{
    cout << "Starting CbmMvdSensorBufferTask::CbmMvdSensorBufferTask() "<< endl;

    fMode          = iMode;
    fBranchName    = "MvdPoint";
    fDigis         = new TClonesArray("CbmMvdDigi");
    fPixelCharge   = new TClonesArray("CbmMvdPixelCharge");
    fPileupManager = NULL;
    fDeltaManager  = NULL;
    fRandGen.SetSeed(2736);
    fEvent       = 0;
    fTime        = 0.;
    fSigmaX      = 0.0005;
    fSigmaY      = 0.0005;
    fNPileup     = 0;
    fNDeltaElect = 0;
    fBgFileName    = "";
    fDeltaFileName = "";
    fBgBufferSize    = 1000;
    fDeltaBufferSize = 10000;
    fReadoutTime = 0.00005;

    fPoints    = new TRefArray();
    //fFluctuate = new MyG4UniversalFluctuationForSi();

    fEpiTh         = 0.0014;
    fSegmentLength = 0.0001;
    fDiffusionCoefficient = 0.0055; // correspondes to the sigma of the gauss with the largest drift length
    fElectronsPerKeV = 276; //3.62 eV deposited for creation of 1 pair electron/hole
    fWidthOfCluster  = 3.5; // in sigmas
    fPixelSizeX      = 0.0030; // in cm
    fPixelSizeY      = 0.0030;
    fCutOnDeltaRays  = 0.00169720; //MeV
    fChargeThreshold = 150; //charge threshold above which the digis are filled
    fFanoSilicium    = 0.115; //not used for the moment, 02 june 08

    fSegmentDepth        = 0;
    fCurrentTotalCharge  = 0;
    fCurrentParticleMass = 0;
    fCurrentParticleMomentum = 0;
    fPixelScanAccelerator    = 0;
    fLandauRandom=new TRandom3();
    /*fLandauMPV=967.;
    fLandauSigma=208;*/

    fShowDebugHistos = kFALSE;

    fPixelSize = 0.0010;
    fPar0 = 3.42157e+02;
    fPar1 = 7.67981e-01;
    fPar2 = 0;


    fLandauMPV=1013.;
    fLandauSigma=159.2;
    fLandauGain=1.56;
    
    /*fLorentzY0=-21.74022;
    fLorentzXc=0.; // -0.04133
    fLorentzW=1.28456;
    fLorentzA=1008.79877; */

    /*
    //line in cluster
    fLorentzY0=-1.7;
    fLorentzXc=0.; // -0.04133
    fLorentzW=1.107;
    fLorentzA=443.35;
    */

    /*
    //column in cluster
    fLorentzY0=-0.16;
    fLorentzXc=0.; // -0.04133
    fLorentzW=0.905;
    fLorentzA=357.02;
    */

    /*
    //params adjusted from "intermediate" mpv values  for mimosa 17
    fLorentzY0=-1.99;
    fLorentzXc=0.; // -0.04133
    fLorentzW=1.056;
    fLorentzA=416.4;
    */

    /*
    //mimosa18 - data from column
    fLorentzY0=-6;
    fLorentzXc=0.;
    fLorentzW=0.95;
    fLorentzA=430.4;
    */

    //mimosa18 - average column + line
    fLorentzY0=-6.1;
    fLorentzXc=0.;
    fLorentzW=1.03;
    fLorentzA=477.2;



    //fLorentzNorm=0.00013010281679422413;
    fLorentzNorm=1;
    
    //SetPixelSize(18);
 }


// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
CbmMvdSensorBufferTask::~CbmMvdSensorBufferTask() {

    if ( fDigis) {
	fDigis->Delete();
	delete fDigis;
    }

    if ( fPileupManager ) delete fPileupManager;
    if ( fDeltaManager  ) delete fDeltaManager;

}




// -----   Virtual public method Exec   ------------------------------------
void CbmMvdSensorBufferTask::Exec(Option_t* opt) {

fPluginHits->Clear("C");

if(fInputPoints->GetEntriesFast() > 0)
{

cout << endl << "Send Input" << endl;
fDetector->SendInput(fInputPoints);
cout << endl << "Execute Chain" << endl;
fDetector->ExecChain();
cout << endl << "End Chain" << endl;


if(bWriteOutput)
    {
    cout << "Start writing Hits" << endl;  
    fPluginHits->AbsorbObjects(fDetector->GetOutputBuffer());
    fPluginHits->Compress();
    fPluginHits->Print();
    cout << endl;
    }
    
}// end of exec
}
// -------------------------------------------------------------------------

Int_t CbmMvdSensorBufferTask::GetMvdGeometry() {
   
  cout << "-I- " << GetName() << " : Reading MVD geometry..." << endl;
  Int_t iStation =  1;
  Int_t volId    = -1;
  fStationMap.clear();

  do 
    {
    // Volume name according to convention
    TString volName  = Form("MVD-S%i-A", iStation);
    //cout << endl << volName << endl;
    volId = gGeoManager->GetUID(volName);
    //cout << endl << volId << endl;
  
    if (volId > -1) 
	{
	
	// Get shape parameters
	TGeoVolume* volume = gGeoManager->GetVolume(volName.Data());   
	TGeoTube* tube = (TGeoTube*) volume->GetShape();     
	//Double_t rmin = tube->GetRmin();   
	//Double_t rmax = tube->GetRmax();
	//Double_t d    = 2. * tube->GetDz();
	
    
	// Full path to node 
	TString nodeName =  Form("cave_1/Beamtimeosetupoobgnum_0/MVD-S%i-AoPartAss_1/MVD-S%i-A_1", iStation, iStation );

	// Get z position of node
	Bool_t nodeFound = gGeoManager->cd(nodeName.Data());
	
	if ( ! nodeFound ) 
	      {
	      cout << "-E- " << GetName() << "::SetMvdGeometry: Node " << nodeName
		   << " not found in geometry!" << endl;
	      Fatal("SetMvdGeometry", "Node not found");
	      }
	Double_t local[3] = {0., 0., 0.};  // Local centre of volume
	Double_t global[3];                // Global centre of volume
	gGeoManager->LocalToMaster(local, global);
	Double_t z = global[2];
	

	// Check for already existing station with the same ID
	// (Just in case, one never knows...)
	if ( fStationMap.find(iStation) != fStationMap.end() )
	      {
	      cout << "-E- " << GetName() << "::GetMvdGeometry: " 
		   << "Volume ID " << iStation << " already in map!" << endl;
	      Fatal("GetMvdGeometry", "Double station number in TGeoManager!");
	      }

	// Create new CbmMvdSensor and add it to the map
	fStationMap[iStation] = new CbmMvdSensor(volName, new CbmMvdMimosa26AHR, volName, nodeName, iStation, volId,0.);
	//fStationMap[iStation]->Print();
      
	iStation++;
	
	}     // Volume found

    } while ( volId > -1 );
    
  return iStation - 1;
}
// -----------------------------------------------------------------------------




// -----   Virtual private method SetParContainers   -----------------------
void CbmMvdSensorBufferTask::SetParContainers() {
    FairRunSim*    sim  = FairRunSim::Instance();
    FairRuntimeDb* rtdb = sim->GetRuntimeDb();
    fGeoPar  = (CbmMvdGeoPar*)  (rtdb->getContainer("CbmMvdGeoPar"));
}
// -------------------------------------------------------------------------


			      
// -----    Virtual private method Init   ----------------------------------
InitStatus CbmMvdSensorBufferTask::Init() {

  cout << "-I- " << GetName() << ": Initialisation..." << endl;


    // Init DebugHistos in case they are needed

    if (fShowDebugHistos) {
	fRandomGeneratorTestHisto = new TH1F("TestHisto","TestHisto",1000,0,12000);
	fResolutionHistoX=new TH1F ("DigiResolutionX","DigiResolutionX", 1000, -.005,.005);
	fResolutionHistoY=new TH1F ("DigiResolutionY","DigiResolutionY", 1000, -.005,.005);
	fPosXY= new TH2F("DigiPointXY","DigiPointXY", 100,-6,6,100,-6,6);
	fpZ= new TH1F ("DigiPointPZ","DigiPointPZ", 1000, -0.5,0.5);
	fPosXinIOut = new TH1F("DigiZInOut","Digi Z InOut", 1000, -0.04,0.04);
	fAngle = new TH1F ("DigiAngle","DigiAngle", 1000,0,90);
	fSegResolutionHistoX= new TH1F("SegmentResolutionX","SegmentResolutionX",1000, -.005,.005);
	fSegResolutionHistoY= new TH1F("SegmentResolutionY","SegmentResolutionY",1000, -.005,.005);
        fSegResolutionHistoZ= new TH1F("SegmentResolutionZ","SegmentResolutionZ",1000, -.005,.005);
	fTotalChargeHisto=new TH1F("TotalChargeHisto","TotalChargeHisto",1000,0,12000);
	fTotalSegmentChargeHisto=new TH1F("TotalSegmentChargeHisto","TotalSegmentChargeHisto",1000,0,12000);

    }




    cout << endl;
    cout << "---------------------------------------------" << endl;
    cout << "-I- Initialising " << GetName() << " ...." << endl;

    // **********  RootManager
    FairRootManager* ioman = FairRootManager::Instance();
    if ( ! ioman ) 
	{
	cout << "-E- " << GetName() << "::Init: No FairRootManager!" << endl;
	return kFATAL;
	}

    // **********  Get input arrays
    fInputPoints = (TClonesArray*) ioman->GetObject("MvdPoint"); 
   
    // **********  Register output array
    fPluginHits = new TClonesArray("CbmMvdHit", 10000);
    ioman->Register("MvdPluginHits", "MvdHit", fPluginHits, kTRUE);
    
    // *****  Get MVD geometry
    Int_t nStations = GetMvdGeometry();
    if ( ! nStations ) 
	{
	cout << "-W- " << GetName() << "::Init: No MVD stations in geometry!"
	     << endl << "   Task will be inactive!" << endl;
        return kERROR;
	}

    // ********** Check for too many pileup / delta events
    if (fNPileup > 200) 
	{
	cout << "-E- " << GetName() << "::Init:  Pileup of " << fNPileup
	     << " too large; maximum buffer size is 1000 events." << endl;
	return kERROR;
	}
	
    if (fNDeltaElect > 20000) 
	{
	cout << "-E- " << GetName() << "::Init:  Delta Pileup of " << fNDeltaElect
	     << " too large; maximum buffer size is 100000 events." << endl;
	return kERROR;
	}

    // **********  Create pileup manager if necessary
    if (fNPileup >= 1 && !(fPileupManager) && fMode == 0 ) 
	{
	  if ( fBgFileName == "" ) 
	      {
	       cout << "-E- " << GetName() << "::Init: Pileup events needed, but no "
		    << " background file name given! " << endl;
	       return kERROR;
	      }
	fPileupManager = new CbmMvdPileupManager(fBgFileName,
						 fBranchName, fBgBufferSize);
	if(fPileupManager->GetNEvents()< 2* fNPileup) 
	    {
	    cout <<"-E- "<< GetName() << ": The size of your BG-File is insufficient to perform the requested pileup" << endl;
	    cout <<"    You need at least events > 2* fNPileup." << endl;
	    cout <<"    Detected: fPileUp = " << fNPileup << ", events in file " << fPileupManager->GetNEvents() << endl; 
	    Fatal("","");  
	    return kERROR;
	    }
	}

    // **********   Create delta electron manager if required
    if (fNDeltaElect >= 1 && !(fDeltaManager) && fMode == 0 ) 
	{
	  if ( fDeltaFileName == "" ) 
	      {
	      cout << "-E- " << GetName() << "::Init: Pileup events needed, but no "
		   << " background file name given! " << endl;
	      return kERROR;
	      }
	  fDeltaManager = new CbmMvdPileupManager(fDeltaFileName,
						fBranchName, fDeltaBufferSize);
	  if(fDeltaManager->GetNEvents()< 2* fNDeltaElect ) 
	      {
	      cout <<"-E- "<< GetName() << ": The size of your Delta-File is insufficient to perform the requested pileup" << endl;
	      cout <<"    You need at least events > 2* fNDeltaElect." << endl;
	      cout <<"    Detected: fNDeltaElect = " << fNDeltaElect << ", events in file " << fDeltaManager->GetNEvents() << endl; 
	      Fatal("","");
	      return kERROR;
	      }
	   
	}

    // Screen output
    cout << GetName() << " initialised with parameters: " << endl;
    //PrintParameters();
    cout << "---------------------------------------------" << endl;

    fDetector = CbmMvdDetector::Instance();
    
    CbmMvdSensorFrameBuffer* frameBuffer = new CbmMvdSensorFrameBuffer;
    CbmMvdPointToHitConverter* converter = new CbmMvdPointToHitConverter;
    CbmMvdSensorTrackingBuffer* trackingBuffer = new CbmMvdSensorTrackingBuffer;
    fDetector->AddPlugin(frameBuffer);
    fDetector->AddPlugin(converter);
    fDetector->AddPlugin(trackingBuffer);
    fDetector->Init();
    
    
    return kSUCCESS;

 

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



// -----   Virtual public method Reinit   ----------------------------------
InitStatus CbmMvdSensorBufferTask::ReInit() {

  // **********  Get MVD geometry
    Int_t nStations = GetMvdGeometry();
    if ( ! nStations ) {
	cout << "-W- " << GetName() << "::ReInit: No MVD stations in geometry!"
	    << endl << "   Task will be inactive!" << endl;
	return kERROR;
    }
    return kSUCCESS;
}
// -------------------------------------------------------------------------



// -----   Virtual method Finish   -----------------------------------------
void CbmMvdSensorBufferTask::Finish() {
/*    cout.setf(ios_base::fixed, ios_base::floatfield);
    cout << endl << "---------------------------------------------" << endl;
    cout << "CbmMvdSensorBufferTask:: Parameters used for digitisation: " << endl;
    cout << "Pixel Size X               : " << setw(8) << setprecision(2)
	<< fPixelSizeX * 10000. << " mum" << endl;
    cout << "Pixel Size Y               : " << setw(8) << setprecision(2)
	<< fPixelSizeY * 10000. << " mum" << endl;
    cout << "Epitaxial layer thickness  : " << setw(8) << setprecision(2)
	<< fEpiTh * 10000. << " mum" << endl;
    cout << "Segment Length             : " << setw(8) << setprecision(2)
	<< fSegmentLength * 10000. << " mum" << endl;
    cout << "Diffusion Coefficient      : " << setw(8) << setprecision(2)
	<< fDiffusionCoefficient * 10000. << " mum" << endl;
    cout << "Width of Cluster           : " << setw(8) << setprecision(2)
	<< fWidthOfCluster << " * sigma " << endl;
    cout << "ElectronsPerKeV 3.62 eV/eh : " << setw(8) << setprecision(2)
	<< fElectronsPerKeV  <<  endl;
    cout << "CutOnDeltaRays             : " << setw(8) << setprecision(8)
	<< fCutOnDeltaRays  << " MeV " <<  endl;
    cout << "ChargeThreshold            : " << setw(8) << setprecision(2)
	<< fChargeThreshold  <<  endl;
    cout << "YOU USED THE LORENTZ DIGITIZER!!! " <<  endl;
    cout << "---------------------------------------------" << endl;
*/


    PrintParameters();


    if (fShowDebugHistos){
	TCanvas* c=new TCanvas("DigiCanvas","DigiCanvas", 150,10,800,600);
	c->Divide(3,3);
	c->cd(1);
        fResolutionHistoX->Draw();
        fResolutionHistoX->Write();
        c->cd(2);
	fResolutionHistoY->Draw();
	fResolutionHistoY->Write();
	c->cd(3);
	fPosXY->Draw();
	fPosXY->Write();
        c->cd(4);
	fpZ->Draw();
	fpZ->Write();
        c->cd(5);
	fPosXinIOut->Draw();
	fPosXinIOut->Write();
        c->cd(6);
	fAngle->Draw();
	fAngle->Write();
	c->cd(7);
	//fSegResolutionHistoX->Draw();
	fSegResolutionHistoX->Write();
	fTotalSegmentChargeHisto->Draw();
	fTotalSegmentChargeHisto->Write();
        c->cd(8);
	fRandomGeneratorTestHisto->Draw();
	fRandomGeneratorTestHisto->Write();

	fSegResolutionHistoY->Write();
        c->cd(9);
	fTotalChargeHisto->Draw();
	fTotalChargeHisto->Write();
	cout << "-I- CbmMvdDigitizerL::Finish - Fit of the total cluster charge"<< endl;
	fTotalChargeHisto->Fit("landau");
	cout << "=============================================================="<< endl; 	 	
	// new TCanvas();
	//fTotalChargeHisto->Draw();
    };

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



// -----   Private method Reset   ------------------------------------------
void CbmMvdSensorBufferTask::Reset() {
    fPluginHits->Delete();

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



// -----   Private method PrintParameters   --------------------------------
void CbmMvdSensorBufferTask::PrintParameters() {
    
    cout.setf(ios_base::fixed, ios_base::floatfield);
    cout << "============================================================" << endl;
    cout << "============== Parameters Buffer Task ======================" << endl;
    cout << "============================================================" << endl;
    
   /** 
    cout << "Pixel Size X               : " << setw(8) << setprecision(2)
	<< fPixelSizeX * 10000. << " mum" << endl;
    cout << "Pixel Size Y               : " << setw(8) << setprecision(2)
	<< fPixelSizeY * 10000. << " mum" << endl;
    cout << "Epitaxial layer thickness  : " << setw(8) << setprecision(2)
	<< fEpiTh * 10000. << " mum" << endl;
    cout << "Segment Length             : " << setw(8) << setprecision(2)
	<< fSegmentLength * 10000. << " mum" << endl;
    cout << "Diffusion Coefficient      : " << setw(8) << setprecision(2)
	<< fDiffusionCoefficient * 10000. << " mum" << endl;
    cout << "Width of Cluster           : " << setw(8) << setprecision(2)
	<< fWidthOfCluster << " * sigma "<< endl;
    cout << "ElectronsPerKeV 3.62 eV/eh : " << setw(8) << setprecision(2)
	<< fElectronsPerKeV  <<  endl;
    cout << "CutOnDeltaRays             : " << setw(8) << setprecision(8)
	<< fCutOnDeltaRays  << " MeV " <<  endl;    cout << "ChargeThreshold            : " << setw(8) << setprecision(2)
	<< fChargeThreshold  <<  endl;
    cout << "Pileup: " << fNPileup << endl;
    cout << "Delta - Pileup: " << fNDeltaElect << endl;**/
    cout << "=============== End Task ===================================" << endl;
 
}
// -------------------------------------------------------------------------  



ClassImp(CbmMvdSensorBufferTask);