//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Uses TPC/CDC tracks with matching Helitron track and determines
//      the difference between the TPC/CDC track extrapolation into the Helitron
//      and the actually measured Helitron hits
//
//
// Environment:
//      Processing of data from FOPI experiment S339
//
// Author List:
//      Paul Buehler                        SMI/OEAW
//
//-----------------------------------------------------------



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

// C++ headers

// Collaborating Class Headers --------
#include "FairRun.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "TVectorD.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TLatex.h"
#include "TLegend.h"

#include "RKTrackRep.h"
#include "GFTrack.h"
#include "GFTrackCand.h"
#include "RecoHits/GFAbsRecoHit.h"
#include "PseudoSpacePoint.h"
#include "GFKalman.h"

#include "FOPIField.h"
#include "FopiEvent.h"
#include "CdcHit.h" 
#include "HeliHit.h" 
#include "HeliTrack.h" 
#include "PlawaTrack.h"
#include "FopiForwardHit.h"

#include "GFException.h"

#define DEBUG 0


//-----------------------------------------------------------
HeliResTask::HeliResTask()
  : fFopiEvBranchName("FopiEvent"),
    fTpcCdcTrackBranchName("TpcCdcPostFit"),
    fHelHitBranchName("HeliHit"),
    fHelTrackBranchName("HeliTrack"),
    fPlaTrackBranchName("PlawaTrack"),
    fTpcHelMatchParsBranchName("TpcCdcHelMatchPars"),
    fHeliResBranchName("HeliResiduals"),
    fminHPmatch(0.),
    fPersistence(kTRUE),
    fVerbose(kFALSE)
{	
  fRad = 3.1415926/180.;

}

//-----------------------------------------------------------
HeliResTask::~HeliResTask(){
}


//-----------------------------------------------------------
InitStatus
HeliResTask::Init() {

  //Get ROOT Manager ------------------
  FairRootManager* ioman= FairRootManager::Instance();
  if(ioman==0)
    {
      Error("HeliResTask::Init","RootManager not instantiated!");
      return kERROR;
    }
  
  // Get input ------------------------
  fFopiEventArray=(TClonesArray*) ioman->GetObject(fFopiEvBranchName);
  if(fFopiEventArray==0) {
    Error("HeliResTask::Init","FOPI event not found!");
    return kERROR;
  } 
  fTpcCdcTracks=(TClonesArray*) ioman->GetObject(fTpcCdcTrackBranchName);
  if(fTpcCdcTracks==0) {
    Error("HeliResTask::Init","TpcCdcTrack array not found!!");
    return kERROR;
  }  
  fHelHitArray=(TClonesArray*) ioman->GetObject(fHelHitBranchName);
  if(fHelHitArray==0) {
    Error("TpcHelMatchingTask::Init","Helitron Hits array not found!");
    return kERROR;
  }  
  fHelTrackArray=(TClonesArray*) ioman->GetObject(fHelTrackBranchName);
  if(fHelTrackArray==0) {
    Error("TpcHelMatchingTask::Init","Helitron track array not found!");
    return kERROR;
  }  
  fPlaTrackArray=(TClonesArray*) ioman->GetObject(fPlaTrackBranchName);
  if(fPlaTrackArray==0) {
    Error("TpcHelMatchingTask::Init","Plawa track array not found!");
    return kERROR;
  }  
  fMatchingPars=(TClonesArray*) ioman->GetObject(fTpcHelMatchParsBranchName);
  if(fMatchingPars==0) {
    Error("TpcHelTrackFittingTask::Init",
      "Matching parameters array not found!!");
    return kERROR;
  }  
  MatchingParsVec = new TVectorD(14);

  // prepare tree for output
  fHeliRes = new TClonesArray("TVector3");
  ioman->Register(fHeliResBranchName,"Tpc",fHeliRes,fPersistence);
  
  counter = 0;

  return kSUCCESS;
}

//-----------------------------------------------------------
void HeliResTask::SetParContainers() {

  //  std::cout<<"TpcResidualTask::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");

}

//-----------------------------------------------------------
void HeliResTask::Exec(Option_t* opt) {
  if (fVerbose) fprintf(stderr,"HeliResTask:Exec\n");

	// Initialisations
  TVector3 pos, mom;
  TMatrixDSym cov;  
  GFDetPlane plane1;
  std::vector<unsigned int> HELcandIDs;
  FopiForwardHit lochit;

  // clean up result buffer
  if(fHeliRes==0) 
    Fatal("HeliResTask::Exec()","No fHeliRes");
  fHeliRes->Delete();

  // get FOPI event
  FopiEvent *fopiEv = (FopiEvent*) fFopiEventArray->At(0);
  if (fopiEv == 0) return;
	
  Int_t nrun = fopiEv->GetRunNum();
	Int_t nspill = fopiEv->GetSpillNum();
	Int_t nevent = fopiEv->GetEvtNum();
  if (fVerbose)
    fprintf(stderr,"Event: %i %i %i\n",nrun,nspill,nevent);
  
  // how many tracks are there
  Int_t NumMatchPars = fMatchingPars->GetEntries();
  if (fVerbose)
    fprintf(stderr,"Number of track combinations: %i\n",NumMatchPars);
  
   // loop over matching Pars
  Int_t ntracks = -1;
  for (Int_t ii=0; ii<NumMatchPars; ii++) {
    MatchingParsVec = (TVectorD*) fMatchingPars->At(ii);
    Double_t *MatchingPars = MatchingParsVec->GetMatrixArray();
    
    // is there are corresponding matched track?
    Int_t trkindex = (Int_t)MatchingPars[0];
    if (trkindex >= 0) {
        
      // load Tpc/Cdc and matching Helitron tracks
      GFTrack *trk = (GFTrack*) fTpcCdcTracks->At((Int_t)MatchingPars[1]);
      Double_t Charge = 1.;
      if (trk->getCharge()<0.) Charge = -1.;
      
      HeliTrack *hsingletrack = (HeliTrack*) fHelTrackArray->At((Int_t)MatchingPars[2]);
      lochit.SetHrsec(hsingletrack->GetHrsec());
      
      // get list of Helitron hits belonging to matching Helitron track
      HELcandIDs.clear(); 
      HELcandIDs = hsingletrack->GetHelIdHits();
      Int_t NumHelTrackHits = HELcandIDs.size();
      
      // loop over Helitron hit points
      for (Int_t jj=0; jj<NumHelTrackHits; jj++) {
        HeliHit *hsinglehit = (HeliHit*) fHelHitArray->At(HELcandIDs[jj]);
        lochit.SetPos(0,hsinglehit->GetHhx(),hsinglehit->GetHhy(),
          hsinglehit->GetHhz());
        

        plane1 = GFDetPlane(
          TVector3(0.,0.,hsinglehit->GetHhz()),
          TVector3(1.,0.,0.),
          TVector3(0.,1.,0.));
        try {
          trk->getPosMomCov(plane1,pos,mom,cov);
        } catch (GFException& e) {
          e.what();
        }        
        
        // fill residual vector into output
        TVector3 residual = TVector3(
          lochit.GetPos(1).X()-pos.X(),
          lochit.GetPos(1).Y()-pos.Y(), Charge*lochit.GetPos(1).Z());
        Int_t entry = fHeliRes->GetEntries();
        new((*fHeliRes)[entry]) TVector3(residual);       
        
        if (fVerbose)
          fprintf(stderr,"Res[%i]: %f %f %f\n",
            jj,lochit.GetPos(1).X()-pos.X(),
            lochit.GetPos(1).Y()-pos.Y(),Charge*lochit.GetPos(1).Z());
      }
      
      // check for PLAWA hit
      if (hsingletrack->GetPHmatch()>fminHPmatch) {
        PlawaTrack *psingletrack = (PlawaTrack*) fPlaTrackArray->At((Int_t)MatchingPars[2]);
        
        Float_t phthe = fRad*psingletrack->GetPhthe();
        Float_t phphi = fRad*psingletrack->GetPhphi();
        Float_t phz   = psingletrack->GetPhz();
        Float_t pr    = phz*tan(phthe);
        
        TVector3 ppos = TVector3(pr*cos(phphi),pr*sin(phphi),phz);

        plane1 = GFDetPlane(
          TVector3(0.,0.,phz),
          TVector3(1.,0.,0.),
          TVector3(0.,1.,0.));
        try {
          trk->getPosMomCov(plane1,pos,mom,cov);
        } catch (GFException& e) {
          e.what();
        }
        
        // fill residual vector into output
        TVector3 residual = TVector3(
          ppos.X()-pos.X(),
          ppos.Y()-pos.Y(), Charge*ppos.Z());
        Int_t entry = fHeliRes->GetEntries();
        new((*fHeliRes)[entry]) TVector3(residual);
        
        if (fVerbose)
          fprintf(stderr,"Res[%i]: %f %f %f\n",
            NumHelTrackHits,ppos.X()-pos.X(),
            ppos.Y()-pos.Y(), Charge*ppos.Z());
      }        
        
    }
  }
  
}

ClassImp(HeliResTask)

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