//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Implementation of class TpcDriftTask
//      see TpcDriftTask.hh for details
//
// Environment:
//      Software developed for the PANDA Detector at FAIR.
//
// Author List:
//      Sebastian Neubert    TUM            (original author)
//
//
//-----------------------------------------------------------

// Panda Headers ----------------------

// This Class' Header ------------------
#include "TpcDriftTask.h"

// C/C++ Headers ----------------------


// Collaborating Class Headers --------
#include "CbmRootManager.h"
#include "TClonesArray.h"
#include "TpcGas.h"
#include "TRandom.h"
#include "TpcPrimaryCluster.h"
#include "TpcDriftedElectron.h"
#include "TpcDevmapCyl.h"
#include "TVector3.h"

// Class Member definitions -----------


TpcDriftTask::TpcDriftTask()
  : CbmTask("TPC Drift"), _persistence(kFALSE), 
    _attach(kTRUE), _diffuseL(kTRUE), _diffuseT(kTRUE), _distort(kTRUE)
{
  _primBranchName = "TpcPrimaryCluster";
  _devFile = "DevMap_29-06-07_E_and_B_new_fieldclass.dat"; //default
  _gas= new TpcGas("NEON-90_CO2-10_B2_PRES1013.asc",400);
  //std::cout<<*_gas<<std::endl;
}


TpcDriftTask::~TpcDriftTask()
{
delete _gas;
delete _devmap;
}

InitStatus
TpcDriftTask::Init()
{
  //Instantiate deviation map
  _devmap = new TpcDevmapCyl(_devFile,_gas->VDrift());
  
  //Get ROOT Manager
  CbmRootManager* ioman= CbmRootManager::Instance();

  if(ioman==0)
    {
      Error("TpcDriftTask::Init","RootManager not instantiated!");
      return kERROR;
    }
  
  // Get input collection
  _primArray=(TClonesArray*) ioman->GetObject(_primBranchName);
  
  if(_primArray==0)
    {
      Error("TpcDriftTask::Init","PrimaryElectron-array not found!");
      return kERROR;
    }
  
  // create and register output array
  _driftedArray = new TClonesArray("TpcDriftedElectron"); 
  ioman->Register("TpcDriftedElectron","Tpc",_driftedArray,_persistence);

  // create histograms with x- and y-shifts
  // these will later be stored in a separate root file 
  _xVariation = new TH1D("xShifts",
			 "x-coordinate shifts from drifting", 200, -0.7, 0.7);
  _yVariation = new TH1D("yShifts",
			 "y-coordinate shifts from drifting" , 200, -0.7, 0.7);
  _xVarAndDriftL = new TH2D("xDrift_vs_DriftLength", 
			    "x-shifts vs. DriftLength",
			    200,-0.7,0.7,200,0,100);
  _yVarAndDriftL = new TH2D("yDrift_vs_DriftLength", 
			    "y-shift vs. DriftLength",
			    200,-0.7,0.7,200,0,100);

  // TODO: Get from geom!!!!
  _zGem=-40;
  

  return kSUCCESS;
}


void
TpcDriftTask::Exec(Option_t* opt)
{
  // Reset output Array
  if(_driftedArray==0) Fatal("TpcPrimCluster::Exec)","No DriftedElectronArray");
  _driftedArray->Delete();

  //loop over incoming electrons
  Int_t nc=_primArray->GetEntriesFast();
  for(int ic=0;ic<nc;++ic){
    TpcPrimaryCluster* pcl=(TpcPrimaryCluster*)_primArray->At(ic);
    //create single electrons
    Int_t q=pcl->q();
    for(Int_t ie=0;ie<q;++ie){
      
      //calculate drift time
      double driftl=pcl->z()-_zGem;
      if(driftl<0)continue;
      //attachment
      if(_attach){
	if ( exp( -driftl * _gas->k() ) < gRandom->Uniform())
	  continue;
      }
      //diffusion
      double dx=0;double dy=0; double dt=0;
      dt=driftl/_gas->VDrift();
      if(_diffuseL){
	double sigmal = _gas->Dl() * sqrt(driftl);
	dt+=gRandom->Gaus(0,sigmal)/_gas->VDrift();
      }
      if(_diffuseT){
	double sigmat = _gas->Dt() * sqrt(driftl);
	dx+=gRandom->Gaus(0,sigmat);
	dy+=gRandom->Gaus(0,sigmat);
      }
      //drift distortions
      if(_distort){
	double posX = pcl->x();
	double posY = pcl->y();
	double posZ = pcl->z();
	TVector3 value = _devmap->value(TVector3(posX, posY, posZ));
	dx+=value.X();
	dy+=value.Y();
	dt+=value.Z() / _gas->VDrift();      
      }
      Int_t size = _driftedArray->GetEntriesFast();
      new((*_driftedArray)[size]) TpcDriftedElectron(pcl->x()+dx,
						       pcl->y()+dy,
						       pcl->t()+dt,
						       pcl);
      //feeding the tracking Histograms with this electrons' data
      FillHistograms(dx, dy, driftl);

    } // end loop over electrons

  } // end loop over clusters 
  std::cout<<_driftedArray->GetEntriesFast()<<" electrons arriving at readout"
	   <<std::endl;
  return;
}

void TpcDriftTask::FillHistograms(double x, double y, double dl) {

  _xVariation->Fill(x);
  _yVariation->Fill(y);
  _xVarAndDriftL->Fill(x, dl);
  _yVarAndDriftL->Fill(y, dl);
}

//WriteHistograms() has to be called once in the runDigi.C macro!
void TpcDriftTask::WriteHistograms() {

  TFile* file = new TFile("histograms.root","RECREATE",
			  "collected data from the simulation process");

  file->mkdir("TpcDriftTask");
  file->cd("TpcDriftTask");

  _xVariation->Write();
  delete _xVariation;
  _yVariation->Write();
  delete _yVariation;
  _xVarAndDriftL->Write();
  delete _xVarAndDriftL;
  _yVarAndDriftL->Write();
  delete _yVarAndDriftL;

  file->Close();
  delete file;
  }

  
ClassImp(TpcDriftTask)