//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Artificially create stuff to test ref-tracking
//      
//
// Environment:
//      Software developed not for PANDA.
//
// Author List:
//      Felix Boehmer                      (original author)
//      Physics Department E18, TUM
//
//
//-----------------------------------------------------------


#include "TpcMickeyMouseMCTask.h"
#include "AbsTpcMMDistortion.h"

#include "TpcCluster.h"
#include "CdcTrack.h"

#include "TRandom.h"
#include "TMath.h"
#include "TVector3.h"
#include <TClonesArray.h>

ClassImp(TpcMickeyMouseMCTask)

TpcMickeyMouseMCTask::TpcMickeyMouseMCTask() 
 :
  fPersistence(kFALSE),
  fClusterOutBranch("TpcCluster"),
  fTrackOutBranch("CdcTrack"),
  fClusterOut(NULL),
  fTrackOut(NULL),
  fSeed(12345),
  fRMin(5.),
  fRMax(15.),
  fMaxTracks(1),
  fMeanClDist(1.),
  fSigRandX(0.),
  fSigRandY(0.),
  fSigRandZ(0.),
  fGlobRot(0)
{
  ;
}

TpcMickeyMouseMCTask::~TpcMickeyMouseMCTask()
{
  ;
}

//almost like your life //
//almost like your endless fight //
//curse the day is long //
//realize you don't belong //
//disconnect somehow //
//never stop complaining now //
//almost like your fight//
//and there it went //
//almost like your life .... //
//that was just your life

InitStatus
TpcMickeyMouseMCTask::Init() {
  
  //connect to the framework:
  FairRootManager* ioman = FairRootManager::Instance();
  if(ioman==0) {
    Fatal("TpcMickeyMouseMCTask::Init()","RootManager not instantiated!");
    return kFATAL;
  }
    
  //register output arrays
  fClusterOut = new TClonesArray("TpcCluster");
  ioman->Register(fClusterOutBranch, "Tpc", fClusterOut, fPersistence);
  //ioman->Register(fClusterOutBranch, "Tpc", fClusterOut, kFALSE);
  fTrackOut = new TClonesArray("CdcTrack");
  //ioman->Register(fTrackOutBranch, "Tpc", fTrackOut, kFALSE);
  ioman->Register(fTrackOutBranch, "Tpc", fTrackOut, fPersistence);
  
  fRand = new TRandom(fSeed);
  
}

void
TpcMickeyMouseMCTask::Exec(Option_t* opt) {

  fClusterOut->Delete();
  fTrackOut->Delete();
  for(unsigned int it=0; it<fTracks.size(); it++) 
    delete fTracks[it];
  fTracks.clear();
  for(unsigned int it=0; it<fClusters.size(); it++) 
    delete fClusters[it];
  fClusters.clear();
  
  //randomize circle(s) intersecting (0,0) 
  unsigned int nTracks;
  if(fMaxTracks==1)
    nTracks=1;
  else
    nTracks = floor(fRand->Uniform()*fMaxTracks) + 1;
 
  //produce CdcTrack object(s)
  for(unsigned int it=0; it<nTracks; it++) {
    double angle = fRand->Uniform()*2*TMath::Pi();
    double rad  = (fRMax-fRMin)*fRand->Uniform();
    TVector3 m;
    m.SetXYZ(1.,0.,0.);
    m.SetPhi(angle);
    m.SetMag(rad);
    fTracks.push_back(new CdcTrack());
    fTracks.back()->SetRadius(rad);
    fTracks.back()->SetMx(m.X());
    fTracks.back()->SetMy(m.Y());
    fTracks.back()->SetValid(1);
    unsigned int itca = fTrackOut->GetEntriesFast();
    CdcTrack* meh = new ((*fTrackOut)[itca]) CdcTrack(*fTracks.back());
  }
  
  unsigned int index = 0;

  //construct tpcClusters along circles 
  for(unsigned int itr=0; itr<fTracks.size(); itr++) {
    CdcTrack* track = fTracks[itr];
    //calculate circumference 
    double L = fTracks[itr]->GetRadius()*2*TMath::Pi();
    //mean distance between clusters:
    int nDiv = floor(L/fMeanClDist);
    
    //get a line from circle center to inner drift cage wall and    
    //step along... this is awkward and quite annoying
    //no better idea right now
    
    //for now: brute force 
    TVector3 r;
    r.SetXYZ(track->GetMx(), track->GetMy(), 0.);
    TVector3 negr=(-1.)*r;
    double theta = r.Theta();
    double dTheta = 0.;
    
            
    //now step along the circle and distribute clusters (Poissonian)
    while(dTheta<2*TMath::Pi()) {
      double step = fRand->Exp(fMeanClDist);
      //corresponding angle:
      //double dt = TMath::Pi()*2/(track->GetRadius()*2*TMath::Pi()/step);
      double dt = step/track->GetRadius();
      dTheta+=dt;
      //check if position is inside active volume
      negr.RotateZ(dt);
      TVector3 pos = r+negr;
      if(pos.Perp()<5. || pos.Perp()>15.)  //the geometrical dimensions
	continue;
      fClusters.push_back(new TpcCluster(pos, 1., index));
      index++;
    }
    
    //shift them following some parametrisation
    for(unsigned int icl=0; icl<fClusters.size(); icl++) {
      TVector3 clpos = (fClusters[icl])->pos();
      //smearing
      double dx = fRand->Gaus(0,fSigRandX);
      double dy = fRand->Gaus(0,fSigRandY);
      double dz = fRand->Gaus(0,fSigRandZ);
      TVector3 dRand(dx,dy,dz);
      clpos+=dRand;
      //rotation
      clpos.RotateZ(fGlobRot);
      (fClusters[icl])->SetPos(clpos);
    }
    
    //local distortions
    for(unsigned int id=0; id<fDistorters.size(); id++)
      fDistorters[id]->distort(fClusters);
    
    //write them out
    for(unsigned int icl=0; icl<fClusters.size(); icl++) 
      TpcCluster* cl = new ((*fClusterOut)[icl]) TpcCluster(*fClusters[icl]);
        
  }

}