//-----------------------------------------------------------
// 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 ----------------------
#include <iostream>
#include <exception>

// Collaborating Class Headers --------
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "TpcDigiPar.h"
#include "TClonesArray.h"
#include "TpcGas.h"
#include "TRandom.h"
#include "TpcPrimaryCluster.h"
#include "TpcDriftedElectron.h"
#include "TpcDevmapCyl.h"
#include "PndDetectorList.h"
#include "TVector3.h"
#include "QAPlotCollection.h"

#include <iostream>
#include <cmath>

using std::cout;
using std::exp;
using std::sqrt;

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


TpcDriftTask::TpcDriftTask()
  : FairTask("TPC Drift", 0), fpersistence(kFALSE),
    fattach(kTRUE), fdiffuseL(kTRUE), fdiffuseT(kTRUE), fdistort(kFALSE), 
    fscale(1), fphicut(kFALSE), finitialized(kFALSE),
    fqa(NULL),fForceManDrift(false)
{
  fprimBranchName = "TpcPrimaryCluster";
  //TODO: parameter management!!!!
  fdevFile = "DevMap_29-06-07_E_and_B_new_fieldclass.dat"; //default
}


TpcDriftTask::~TpcDriftTask()
{
if(fdevmap!=NULL)delete fdevmap;
}

InitStatus
TpcDriftTask::Init()
{ 

  finitialized=false;

  //Get ROOT Manager
  FairRootManager* ioman= FairRootManager::Instance();

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

  // create histograms with x- and y-shifts
  // these will later be stored in a separate root file
  if(fqa!=NULL){
  fxVariation = fqa->getTH1D("xShifts",
			 "x-coordinate shifts from drifting", 200, -0.7, 0.7);
  fyVariation = fqa->getTH1D("yShifts",
			 "y-coordinate shifts from drifting" , 200, -0.7, 0.7);
  fxVarAndDriftL = fqa->getTH2D("xDrift_vs_DriftLength", 
			    "x-shifts vs. DriftLength",
			    200,-0.7,0.7,200,0,100);
  fyVarAndDriftL = fqa->getTH2D("yDrift_vs_DriftLength", 
			    "y-shift vs. DriftLength",
			    200,-0.7,0.7,200,0,100);
  }

  fgas=fpar->getGas();
    fDriftVel = fpar->getManDriftVel();
  if(fDriftVel>0){
    fForceManDrift=true;
  }  else{
    fDriftVel = fgas->VDrift();
    fForceManDrift=false;
  }
      
  fdiffuseL=fpar->getDiffuseL();
  fdiffuseT=fpar->getDiffuseT();
  fattach=fpar->getAttach();
  
//Instantiate deviation map
  if(fdistort){
    fdevmap = new TpcDevmapCyl(fdevFile,fDriftVel);
    fdevmap->setScale(fscale);
    if(!fdevmap->loaded()) {
      Fatal("TpcDriftTask::Init","Deviation Map not loaded! Aborting...");
      return kERROR;
    }
  }
  
  finitialized=true;
  return kSUCCESS;
}

void 
TpcDriftTask::SetParContainers() {

  std::cout<<"TpcDriftTask::SetParContainers"<<std::endl;
  std::cout.flush();

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

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

  // Get Tpc digitisation parameter container
  fpar= (TpcDigiPar*) db->getContainer("TpcDigiPar");
  if (! fpar ) Fatal("SetParContainers", "TpcDigiPar not found");
}



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

  //loop over incoming electrons
  Int_t nc=fprimArray->GetEntriesFast();

  std::vector<TVector3>diffVec;
  //if(fdiffuseL and fdiffuseT)
	  //diffVec=MakeDiffusion(nc);

  for(int ic=0;ic<nc;++ic){
    TpcPrimaryCluster* pcl=(TpcPrimaryCluster*)fprimArray->At(ic);

    if(fphicut){
      double phi=pcl->pos().Phi();
      if(phi<fphimin || phi>fphimax)continue;
    }
    //create single electrons
    Int_t q=pcl->q();
    for(Int_t ie=0;ie<q;++ie){
      
      //calculate drift time
      double driftl=pcl->z();
      if(driftl<0)continue;
      //attachment
      if(fattach){
        if ( exp( -driftl * fgas->k() ) < gRandom->Uniform())
          continue;
      }
      //diffusion
      double dx=0;double dy=0; double dt=0; double t=0;
      t=driftl/fDriftVel;
      if(fdiffuseL)// and not fdiffuseT)
    	  {
        double sigmal = fgas->Dl() * sqrt(driftl);
        dt+=gRandom->Gaus(0,sigmal)/fDriftVel;
      }
      if(fdiffuseT)// and not fdiffuseL)
    	  {
        double sigmat = fgas->Dt() * sqrt(driftl);
        dx+=gRandom->Gaus(0,sigmat);
        dy+=gRandom->Gaus(0,sigmat);
      }
      //drift distortions
      if(fdistort){
        double posX = pcl->x();
        double posY = pcl->y();
        double posZ = pcl->z();
        TVector3 value = fdevmap->value(TVector3(posX, posY, posZ));
        dx+=value.X();
        dy+=value.Y();
        dt+=value.Z() / fDriftVel;
      }
      Int_t size = fdriftedArray->GetEntriesFast();
      TpcDriftedElectron* myElectron = new((*fdriftedArray)[size]) TpcDriftedElectron(pcl->x()+dx,
                                                                                      pcl->y()+dy,
                                                                                      pcl->t()+dt+t,
                                                                                      pcl);
      myElectron->setIndex(size);
      myElectron->setDist(dx,dy,dt);

      //feeding the tracking Histograms with this electrons' data
      FillHistograms(dx, dy, driftl);
      //std::cout<<"x="<<pcl->x()<<"  y="<< pcl->y()<<std::endl;
    } // end loop over electrons

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

void TpcDriftTask::FillHistograms(double x, double y, double dl) {
	if(fqa!=NULL){
		fxVariation->Fill(x);
		fyVariation->Fill(y);
		fxVarAndDriftL->Fill(x, dl);
		fyVarAndDriftL->Fill(y, dl);
	}
}

//WriteHistograms() has to be called once in the runDigi.C macro!
void TpcDriftTask::WriteHistograms() {
  if(!finitialized)return;
  TFile* file=FairRootManager::Instance()->GetOutFile();

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

  fxVariation->Write();
  delete fxVariation;
  fyVariation->Write();
  delete fyVariation;
  fxVarAndDriftL->Write();
  delete fxVarAndDriftL;
  fyVarAndDriftL->Write();
  delete fyVarAndDriftL;

 }

  
ClassImp(TpcDriftTask)