#include "PndLheHitsMaker.h"

#include "PndLheHit.h"
#include "PndLheCandidate.h"

#include "FairRunAna.h"
#include "PndTpcPoint.h"
#include "PndTpcCluster.h"
#include "PndSttPoint.h"
#include "PndSttHit.h"
#include "PndSttHelixHit.h"
#include "PndSdsMCPoint.h"
#include "PndSdsHit.h"
#include "PndSdsCluster.h"
#include "PndSdsDigiPixel.h"
#include "PndSdsDigiStrip.h"
#include "PndEmcCluster.h"
#include "PndEmcBump.h"
#include "PndGemMCPoint.h"
#include "PndGemHit.h"
#include "PndMCTrack.h"
#include "FairMCApplication.h"
#include "FairRootManager.h"
#include "FairTask.h"
#include "TRandom.h"

#include <iomanip>

//_________________________________________________________________
PndLheHitsMaker::PndLheHitsMaker() {
  //---

  fLheHits  = new TClonesArray("PndLheHit");
  fGeantTracks = new TClonesArray("PndLheCandidate");
  fMvdMode       = 1;   
  fMvdSimMode    = 1; 
  fTpcMode       = 1; 
  fSttMode       = 0;
  fSttSimMode    = 1; 
  fEmcMode       = 0;
  fGemMode       = 0;
  fGemSimMode    = 1;  
  fSimulation    = kTRUE;
  fPersistence   = kTRUE;
  fTpcResolution = -1.;
  fMvdResolution = -1.;
  fSttResolution = -1.;
  fGemResolution = -1.;
  fVerbose       = kFALSE;
}

//_________________________________________________________________
PndLheHitsMaker::PndLheHitsMaker(const char *name,
				       const char *title):FairTask(name) {
  //---
  fLheHits = new TClonesArray("PndLheHit");
  fGeantTracks = new TClonesArray("PndLheCandidate");
  fMvdMode       = 1; 
  fMvdSimMode    = 1; 
  fTpcMode       = 1;
  fSttMode       = 0;
  fSttSimMode    = 1; 
  fEmcMode       = 0; 
  fGemMode       = 0; 
  fGemSimMode    = 1;  
  fSimulation    = kTRUE;
  fPersistence   = kTRUE;
  fTpcResolution = -1.;
  fMvdResolution = -1.;
  fSttResolution = -1.;
  fGemResolution = -1.;
  fVerbose       = kFALSE;
}

//_________________________________________________________________
PndLheHitsMaker::~PndLheHitsMaker() {
  //
  FairRootManager *fManager =FairRootManager::Instance();
  fManager->Write();
}

//_________________________________________________________________
InitStatus PndLheHitsMaker::ReInit() {
  // --- Update the pointers after reinitialisation

  return kSUCCESS;
}

//_________________________________________________________________
InitStatus PndLheHitsMaker::Init() {

  //  cout << "InitStatus PndLheHitsMaker::Init\n\n";
  cout << "-I- PndLheHitsMaker::Init" << endl;
  
  fCuts = PndLheTrackCuts::Instance();

  FairRootManager *fManager =FairRootManager::Instance();	
  
  switch (fTpcMode) 
    {
    case 0:
      cout << "-I- PndLheHitsMaker::Init: No TPC detector is used" << endl;
      break;
      
    case 1:
      fTpcInput   = (TClonesArray *)fManager->GetObject("PndTpcPoint");
      if ( ! fTpcInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No TpcPoint array! Switching TPC OFF" << endl;
	  fTpcMode = 0;
	}
      else
	{
	  if (fTpcResolution<0.) cout << "-I- PndLheHitsMaker::Init: Using PndTpcPoint, no position smearing" << endl;
	  else cout << "-I- PndLheHitsMaker::Init: Using PndTpcPoint, position smearing " << fTpcResolution << " [cm]" << endl;
	}
      break;
      
    case 2:
      fTpcInput   = (TClonesArray *)fManager->GetObject("PndTpcCluster");
      if ( ! fTpcInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No TpcCluster array! Switching TPC OFF" << endl;
	  fTpcMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using PndTpcCluster" << endl;
	}
      break;

    default:
      cout << "-E- PndLheHitsMaker::Init: Wrong TPC mode. Switching TPC OFF" << endl;
      fTpcMode = 0;
    }
  
 switch (fSttMode) 
    {
    case 0:
      cout << "-I- PndLheHitsMaker::Init: No STT detector is used" << endl;
      break;
      
    case 1:
      fSttMCArray   = (TClonesArray *)fManager->GetObject("STTPoint");
      if ( ! fSttMCArray ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No SttPoint array! Switching STT OFF" << endl;
	  fSttMode = 0;
	}
      else
	{
	  if (fSttResolution<0.) cout << "-I- PndLheHitsMaker::Init: Using PndSttPoint, no position smearing" << endl;
	  else cout << "-I- PndLheHitsMaker::Init: Using PndSttPoint, position smearing " << fSttResolution << " [cm]" << endl;
	}
      break;
      
    case 2:
      fSttHitInput   = (TClonesArray *)fManager->GetObject("STTHit");
      if ( ! fSttHitInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No STTHit array! Switching STT OFF" << endl;
	  fSttMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using PndSttHit" << endl;
	}
      break;
      
    case 3:
      fSttHelixInput   = (TClonesArray *)fManager->GetObject("SttHelixHit");
      if ( ! fSttHelixInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No SttHelixHit array! Switching STT OFF" << endl;
	  fSttMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using PndSttHelixHit" << endl;
	} 

      fSttMCArray = (TClonesArray*) fManager->GetObject("STTPoint");
      if ( ! fSttMCArray )         fSttSimMode = 0; 
      
      break;
      
    case 4:
      fSttHitInput   = (TClonesArray *)fManager->GetObject("STTHit");
      if ( ! fSttHitInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No STTHit array! Switching STT OFF" << endl;
	  fSttMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using PndSttHit" << endl;
	}
      fSttHelixInput   = (TClonesArray *)fManager->GetObject("SttHelixHit");
      if ( ! fSttHelixInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No SttHelixHit array! Switching STT OFF" << endl;
	  fSttMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using PndSttHelixHit" << endl;
	}
      
      fSttMCArray = (TClonesArray*) fManager->GetObject("STTPoint");
      if ( ! fSttMCArray )         fSttSimMode = 0; 
      
      break;
      
    case 5:
      fSttPRHelixInput   = (TClonesArray *)fManager->GetObject("STTPRHelixHit");
      if ( ! fSttPRHelixInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No STTPRHelixHit array! Switching STT OFF" << endl;
	  fSttMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using PndSttPRHelixHit" << endl;
	}

      fSttMCArray = (TClonesArray*) fManager->GetObject("STTPoint");
      if ( ! fSttMCArray )         fSttSimMode = 0;
      break;
      
    default:
      cout << "-E- PndLheHitsMaker::Init: Wrong STT mode. Switching STT OFF" << endl;
      fSttMode = 0;
    }
 
 if (fSttMode>0 && fTpcMode > 0)
   {
     cout << "-E- PndLheHitsMaker::Init: Both STT and TPC mode. Too many tracking detectors for my tastes!!!" << endl;
     return kFATAL;
   }
 

  switch (fMvdMode) 
    {
    case 0:
      cout << "-I- PndLheHitsMaker::Init: No MVD detector is used" << endl;
      break;

    case 1:
      fMvdMCArray = (TClonesArray*) fManager->GetObject("MVDPoint");
      if ( ! fMvdMCArray ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No MVDPoint array! Switching MVDPoint OFF" << endl;
	  fMvdMode = 0;
	}
      else
	{
	  if (fMvdResolution<0.) cout << "-I- PndLheHitsMaker::Init: Using MVDPoint, no position smearing" << endl;
	  else cout << "-I- PndLheHitsMaker::Init: Using MVDPoint, position smearing " << fMvdResolution << " [cm]" << endl;
	}
      break;
      
    case 2: 
      fMvdStripHitArray = (TClonesArray*) fManager->GetObject("MVDHitsStrip");
      if ( !fMvdStripHitArray)
	{
	  cout << "-W- PndLheHitsMaker::Init: No MVDHitsStrip array!" << endl;
	}
      fMvdPixelHitArray = (TClonesArray*) fManager->GetObject("MVDHitsPixel");
      if ( !fMvdPixelHitArray)
	{
	  cout << "-W- PndMvdIdealTrackFinderTask::Init: " << "No fMvdPixelHitArray!" << endl;
	}
      
      if (( ! fMvdStripHitArray ) &&  ( ! fMvdPixelHitArray ))
	{
	  cout << "-W- PndLheHitsMaker::Init: No MVD hits array! Switching MVD OFF" << endl;
	  fMvdMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using MVDHit" << endl;
	}
     
      fMvdMCArray = (TClonesArray*) fManager->GetObject("MVDPoint");
      if ( ! fMvdMCArray )         fMvdSimMode = 0;

      if ( fMvdSimMode == 0)  cout << "-W- PndLheHitsMaker::Init: No MVD MC information is available" << endl;
      
      break;
      
    default:
      cout << "-E- PndLheHitsMaker::Init: Wrong MVD mode. Switching MVD OFF" << endl;
      fMvdMode = 0;
    }

  switch (fEmcMode) 
    {
    case 0:
      cout << "-I- PndLheHitsMaker::Init: No EMC detector is used" << endl;
      break;
      
    case 2:
      fEmcInput = (TClonesArray*) fManager->GetObject("EmcCluster");
      if ( ! fEmcInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No EmcCluster array! Switching EmcCluster OFF" << endl;
	  fEmcMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using EmcCLuster" << endl;
	}
      break;
      
    case 3:
      fEmcInput = (TClonesArray*) fManager->GetObject("EmcBump");
      if ( ! fEmcInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No EmcBump array! Switching EmcBump OFF" << endl;
	  fEmcMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using EmcBump" << endl;
	}
      break;
      
    default:
      cout << "-E- PndLheHitsMaker::Init: Wrong EMC mode. Switching EMC OFF" << endl;
      fEmcMode = 0;
    }
  
  switch (fGemMode) 
    {
    case 0:
      cout << "-I- PndLheHitsMaker::Init: No GEM detector is used" << endl;
      break;
      
    case 1:
      fGemInput   = (TClonesArray *)fManager->GetObject("GEMPoint");
      if ( ! fGemInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No GEMPoint array! Switching GEM OFF" << endl;
	  fGemMode = 0;
	}
      else
	{
	  if (fGemResolution<0.) cout << "-I- PndLheHitsMaker::Init: Using GEMPoint, no position smearing" << endl;
	  else cout << "-I- PndLheHitsMaker::Init: Using GEMPoint, position smearing " << fGemResolution << " [cm]" << endl;
	}
      break;
      
    case 2:
      fGemInput   = (TClonesArray *)fManager->GetObject("GEMHit");
      if ( ! fGemInput ) 
	{
	  cout << "-W- PndLheHitsMaker::Init: No GEMHit array! Switching GEM OFF" << endl;
	  fGemMode = 0;
	}
      else
	{
	  cout << "-I- PndLheHitsMaker::Init: Using GEMHit" << endl;
	}
      
      fGemMCArray = (TClonesArray*) fManager->GetObject("GEMPoint");
      if ( ! fGemMCArray )         fGemSimMode = 0; 
      break;
      
    default:
      cout << "-E- PndLheHitsMaker::Init: Wrong GEM mode. Switching GEM OFF" << endl;
      fGemMode = 0;
    }
  
  if ((fMvdMode+fTpcMode+fSttMode+fEmcMode+fGemMode)==0)
    {
      cout << "-E- PndLheHitsMaker::Init: No active detectors! Exit" << endl;
      return kFATAL;
    }
  
  if (fSimulation)
    { 
      fListMCtracks = (TClonesArray *)fManager->GetObject("MCTrack");
      if ( ! fListMCtracks ) {
	cout << "-W- PndLheHitsMaker::Init: No MCTrack array!" << "\t" << "Switching Simulation OFF" << endl;
	fSimulation = kFALSE;
      }
      else cout << "-I- PndLheHitsMaker::Init: Using MCTrack" << endl;
    }
  
  if (fVerbose)   cout << "-I- PndLheHitsMaker::Init: Verbose Mode ON" << endl;
  
  fNofEvents = 0;
  
  Register();
  
  return kSUCCESS;
  
}

//_________________________________________________________________
void PndLheHitsMaker::GetMvdPoints() {
  // Taking points from PndSdsMCPoint
  
  if (fVerbose)
    cout << " -I- PndLheHitsMaker::GetMvdPoints(): MVD points entries " << fMvdMCArray->GetEntriesFast() <<endl;
  
  for (int j=0; j < fMvdMCArray->GetEntriesFast(); j++ ) 
    {
      PndSdsMCPoint* point = (PndSdsMCPoint*) fMvdMCArray->At(j);
      
      if (fVerbose) point->Print(" "); //PR(point->GetTrackID());
      
      PndLheHit* hit = AddHit();
      hit->SetHitNumber(fNHit++);
      
      if (fTpcResolution<=0.)
	{
	  hit->SetX(0.5*(point->GetX()+point->GetXOut()));
	  hit->SetY(0.5*(point->GetY()+point->GetYOut()));
	  hit->SetZ(0.5*(point->GetZ()+point->GetZOut()));
	  
	  if (fVerbose) cout << "MVD POINT " 
			     << hit->GetX() << " " 
			     << hit->GetY() << " " 
			     << hit->GetZ() << " RADIUS " 
			     << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";	
	  
	  hit->SetDx(0.005);
	  hit->SetDy(0.005);
	  hit->SetDz(0.005);
	}
      else
	{
	  hit->SetX(gRandom->Gaus(point->GetX(),fMvdResolution));
	  hit->SetY(gRandom->Gaus(point->GetY(),fMvdResolution));
	  hit->SetZ(gRandom->Gaus(point->GetZ(),fMvdResolution));
	  
	  if (fVerbose) cout << "MVD POINT SMEARED " 
			     << hit->GetX() << " " 
			     << hit->GetY() << " " 
			     << hit->GetZ() << " RADIUS " 
			     << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
	  
	  hit->SetDx(fMvdResolution);
	  hit->SetDy(fMvdResolution);
	  hit->SetDz(fMvdResolution);
	}
      
      hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("MVDPoint"));
      hit->SetTrackID(point->GetTrackID());
      hit->SetRefIndex(j);
      hit->SetLink(FairLink("MVDPoint", j));
      
      if (fVerbose)  hit->Print();
    } // end of MVDPoint loop  
}

//_________________________________________________________________
void PndLheHitsMaker::GetMvdHits() {
  // Taking points from PndMvdHits
  
  for (int j2=0; j2 < fMvdPixelHitArray->GetEntriesFast(); j2++ ) 
    {
      PndSdsHit* point = (PndSdsHit*) fMvdPixelHitArray->At(j2);
      
      PndLheHit* hit = AddHit();
      hit->SetHitNumber(fNHit++);
      
      hit->SetX(point->GetX());
      hit->SetY(point->GetY());
      hit->SetZ(point->GetZ());
      
      if (fVerbose) cout << "MVD PIXEL HIT " 
			 << hit->GetX() << " " 
			 << hit->GetY() << " " 
			 << hit->GetZ() << " RADIUS " 
			 << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";	
      
      // hit->SetDx(point->GetDx());
      // hit->SetDy(point->GetDy());
      // hit->SetDz(point->GetDz());
     
      hit->SetDx(0.005);
      hit->SetDy(0.005);
      hit->SetDz(0.005);

      hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("MVDHitsPixel"));
      
      if ( (fMvdSimMode) && (point->GetRefIndex()!=-1) )
	{
	  PndSdsMCPoint* myPoint = (PndSdsMCPoint*)(fMvdMCArray->At(point->GetRefIndex()));
	  hit->SetTrackID(myPoint->GetTrackID());
	}
      hit->SetRefIndex(j2);
      hit->SetLink(FairLink("MVDHitsPixel", j2));
      if (fVerbose)  hit->Print();
    } // end of MVDHitsPixel loop  
  

 
  for (int j=0; j < fMvdStripHitArray->GetEntriesFast(); j++ ) 
    {
      PndSdsHit* point = (PndSdsHit*) fMvdStripHitArray->At(j);
      
      PndLheHit* hit = AddHit();
      hit->SetHitNumber(fNHit++);
      
      hit->SetX(point->GetX());
      hit->SetY(point->GetY());
      hit->SetZ(point->GetZ());
      
      if (fVerbose) cout << "MVD STRIP HIT " 
			 << hit->GetX() << " " 
			 << hit->GetY() << " " 
			 << hit->GetZ() << " RADIUS " 
			 << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";	
      
      // hit->SetDx(point->GetDx());
      // hit->SetDy(point->GetDy());
      // hit->SetDz(point->GetDz());
       
      hit->SetDx(0.005);
      hit->SetDy(0.005);
      hit->SetDz(0.005);
      hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("MVDHitsStrip"));

      if ( (fMvdSimMode) && (point->GetRefIndex()!=-1) )
	{
	  PndSdsMCPoint* myPoint = (PndSdsMCPoint*)(fMvdMCArray->At(point->GetRefIndex()));
	  hit->SetTrackID(myPoint->GetTrackID());
	}
      hit->SetRefIndex(j);
      hit->SetLink(FairLink("MVDHitsStrip", j));
      if (fVerbose)  hit->Print();
    } // end of MVDHitsStrip loop  

}
  
//_________________________________________________________________
void PndLheHitsMaker::GetTpcPoints() {
   // Taking points from PndTpcPoint
  
  if (fVerbose)
    cout << " PndLheHitsMaker::GetTpcHits(): Tpc points entries " << fTpcInput->GetEntriesFast() <<endl;
  
  
  for (int j=0; j < fTpcInput->GetEntriesFast(); j++ ) {
    PndTpcPoint* point = (PndTpcPoint*) fTpcInput->At(j);
    
    if (fVerbose) point->Print(" "); //PR(point->GetTrackID());

    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    if (fTpcResolution<=0.)
      {
	hit->SetX(point->GetX());
	hit->SetY(point->GetY());
	hit->SetZ(point->GetZ());
	
	if (fVerbose) cout << "TPC POINT " 
			   << hit->GetX() << " " 
			   << hit->GetY() << " " 
			   << hit->GetZ() << " RADIUS " 
			   << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
	
	hit->SetDx(0.03);
	hit->SetDy(0.03);
	hit->SetDz(0.03);
      }
    else
      {
	hit->SetX(gRandom->Gaus(point->GetX(),fTpcResolution));
	hit->SetY(gRandom->Gaus(point->GetY(),fTpcResolution));
	hit->SetZ(gRandom->Gaus(point->GetZ(),2*fTpcResolution));
	
	if (fVerbose) cout << "TPC POINT SMEARED " 
			   << hit->GetX() << " " 
			   << hit->GetY() << " " 
			   << hit->GetZ() << " RADIUS " 
			   << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
	
	hit->SetDx(fTpcResolution);
	hit->SetDy(fTpcResolution);
	hit->SetDz(2*fTpcResolution);
      }
    
    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("PndTpcPoint"));
    hit->SetTrackID(point->GetTrackID());
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("PndTpcPoint", j));
    if (fVerbose)  hit->Print();
  }  // end of TpcPoints loop
  
} 

//_________________________________________________________________
void PndLheHitsMaker::GetTpcClusters() {
   // Taking points from PndTpcCluster

  for (int j=0; j < fTpcInput->GetEntriesFast(); j++ ) {
    PndTpcCluster* clu = (PndTpcCluster*) fTpcInput->At(j);
    
    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(clu->pos().X());  
    hit->SetY(clu->pos().Y());
    hit->SetZ(clu->pos().Z());
  
    
    if (fVerbose) cout << "TPC CLU " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(clu->sig().X());  
    hit->SetDy(clu->sig().Y());
    hit->SetDz(clu->sig().Z());

    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("PndTpcCluster"));
    hit->SetTrackID(clu->mcId().DominantID().mctrackID());
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("PndTpcCluster", j));
    if (fVerbose)  hit->Print();
    
  }  // end of TpcCluster loop  
} 
//_________________________________________________________________
void PndLheHitsMaker::GetSttPoints() {
   // Taking points from PndSttPoint
  
  if (fVerbose)
    cout << " PndLheHitsMaker::GetSttPoints(): Stt points entries " << fSttMCArray->GetEntriesFast() <<endl;
  
  
  for (int j=0; j < fSttMCArray->GetEntriesFast(); j++ ) {
    PndSttPoint* point = (PndSttPoint*) fSttMCArray->At(j);
    
    if (fVerbose) point->Print(" "); //PR(point->GetTrackID());

    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    if (fSttResolution<=0.)
      {
	hit->SetX(point->GetX());
	hit->SetY(point->GetY());
	hit->SetZ(point->GetZ());
	
	if (fVerbose) cout << "STT POINT " 
			   << hit->GetX() << " " 
			   << hit->GetY() << " " 
			   << hit->GetZ() << " RADIUS " 
			   << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
	
	hit->SetDx(.5);
	hit->SetDy(.5);
	hit->SetDz(.5);
      }
    else
      {
	hit->SetX(gRandom->Gaus(point->GetX(),fSttResolution));
	hit->SetY(gRandom->Gaus(point->GetY(),fSttResolution));
	hit->SetZ(gRandom->Gaus(point->GetZ(),2*fSttResolution));
	
	if (fVerbose) cout << "STT POINT SMEARED " 
			   << hit->GetX() << " " 
			   << hit->GetY() << " " 
			   << hit->GetZ() << " RADIUS " 
			   << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
	
	hit->SetDx(fSttResolution);
	hit->SetDy(fSttResolution);
	hit->SetDz(2*fSttResolution);
      }
    
    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("STTPoint"));
    hit->SetTrackID(point->GetTrackID());
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("STTPoint", j));
    if (fVerbose)  hit->Print();
  }  // end of SttPoints loop
  
}


//_________________________________________________________________
void PndLheHitsMaker::GetSttHit() {
   // Taking points from PndSttHit

  for (int j=0; j < fSttHitInput->GetEntriesFast(); j++ ) {
    PndSttHit* sttHit = (PndSttHit*) fSttHitInput->At(j);
    
    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(sttHit->GetX());  
    hit->SetY(sttHit->GetY());
    hit->SetZ(sttHit->GetZ());
  
    
    if (fVerbose) cout << "STT HIT " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(.5);  
    hit->SetDy(.5);
    hit->SetDz(.5);

    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("STTHit"));
    hit->SetTrackID(hit->GetRefIndex());
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("STTHit", j));
    
    if (fVerbose)  hit->Print();
    
  }  // end of SttHit loop  
} 

//_________________________________________________________________
void PndLheHitsMaker::GetSttHelixHit() {
   // Taking points from PndSttHelixHit

  for (int j=0; j < fSttHelixInput->GetEntriesFast(); j++ ) {
    PndSttHelixHit* sttHit = (PndSttHelixHit*) fSttHelixInput->At(j);
     if ((isnan(sttHit->GetX())) || (isnan(sttHit->GetZ())))
       {
        sttHit->Print();
        cout << "-W- PndLheHitsMaker::GetSttHelixHit: SttHelixHit with NaN - skipped (requires STT bug fix!!)" << endl;
        continue;
      }
      if (sttHit->GetZ()==0)
      {
        sttHit->Print();
        cout << "-W- PndLheHitsMaker::GetSttHelixHit: SttHelixHit with GetZ() == 0 - skipped (requires STT bug fix!!)" << endl;
        continue;
      }

    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(sttHit->GetX());  
    hit->SetY(sttHit->GetY());
    hit->SetZ(sttHit->GetZ());
  
    
    if (fVerbose) cout << "STT HELIX HIT " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(.5);  
    hit->SetDy(.5);
    hit->SetDz(.5);

    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("STTHelixHit"));
    hit->SetTrackID(-1);
    if ( (fSttSimMode) && (sttHit->GetRefIndex()!=-1) )
    {
		PndSttPoint* myPoint = (PndSttPoint*)(fSttMCArray->At(sttHit->GetRefIndex()));
		hit->SetTrackID(myPoint->GetTrackID());
    }
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("STTHelixHit", j));

    if (fVerbose)  hit->Print();
  }  // end of SttHelixHit loop  
} 

//_________________________________________________________________
void PndLheHitsMaker::GetSttPRHelixHit() {
   // Taking points from PndSttHelixHit

  for (int j=0; j < fSttPRHelixInput->GetEntriesFast(); j++ ) {
    PndSttHelixHit* sttHit = (PndSttHelixHit*) fSttPRHelixInput->At(j);
    if ((isnan(sttHit->GetX())) || (isnan(sttHit->GetZ())))
      {
        sttHit->Print();
        cout << "-W- PndLheHitsMaker::GetSttPRHelixHit: SttHelixHit with NaN - skipped (requires STT bug fix!!)" << endl;
        continue;
      }
      if (sttHit->GetZ()==0)
      {
        sttHit->Print();
        cout << "-W- PndLheHitsMaker::GetSttPRHelixHit: SttHelixHit with GetZ() == 0 - skipped (requires STT bug fix!!)" << endl;
        continue;
      }

    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(sttHit->GetX());  
    hit->SetY(sttHit->GetY());
    hit->SetZ(sttHit->GetZ());
  
    
    if (fVerbose) cout << "STT HELIX HIT " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(.5);  
    hit->SetDy(.5);
    hit->SetDz(.5);

    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("STTHelixHit"));
    hit->SetTrackID(-1);
    if ( (fSttSimMode) && (sttHit->GetRefIndex()!=-1) )
    {
		PndSttPoint* myPoint = (PndSttPoint*)(fSttMCArray->At(sttHit->GetRefIndex()));
		hit->SetTrackID(myPoint->GetTrackID());
    }
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("STTHelixHit", j));

    if (fVerbose)  hit->Print();
  }  // end of SttHelixHit loop  
} 
//_________________________________________________________________
void PndLheHitsMaker::GetSttHelixHitMC() {
   // Taking points from PndSttHelixHit

  for (int j=0; j < fSttHelixInput->GetEntriesFast(); j++ ) {
    PndSttHelixHit* sttHelixHit = (PndSttHelixHit*) fSttHelixInput->At(j);
    PndSttHit* sttHit = (PndSttHit*) fSttHitInput->At(sttHelixHit->GetHitIndex());
    PndSttPoint* sttPoint = (PndSttPoint*) fSttMCArray->At(sttHit->GetRefIndex());
   
    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(sttPoint->GetXtot());  
    hit->SetY(sttPoint->GetYtot());
    hit->SetZ(sttPoint->GetZtot());
  
    
    if (fVerbose) cout << "STT HELIX HIT " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(.5);  
    hit->SetDy(.5);
    hit->SetDz(.5);

    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("STTHelixHit"));
    hit->SetTrackID(sttPoint->GetTrackID());
    
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("STTHelixHit", j));

    if (fVerbose)  hit->Print();
  }  // end of SttHelixHitMC loop  
} 
//_________________________________________________________________
void PndLheHitsMaker::GetEmcClusters() {
   // Taking points from PndEmcCluster

  for (int j=0; j < fEmcInput->GetEntriesFast(); j++ ) {
    PndEmcCluster* clu = (PndEmcCluster*) fEmcInput->At(j);
    
    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(clu->x());  
    hit->SetY(clu->y());
    hit->SetZ(clu->z());
  
    
    if (fVerbose) cout << "EMC CLU " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(1.);  
    hit->SetDy(1.);
    hit->SetDz(1.);

    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("EmcCluster"));
    //hit->SetTrackID(point->GetTrackID());
    hit->SetTrackID(-1);
    hit->SetRefIndex(j);

    hit->SetLink(FairLink("EmcCluster", j));
    if (fVerbose)  hit->Print();
    
  }  // end of EmcCluster loop  
} 
//_________________________________________________________________
void PndLheHitsMaker::GetEmcBumps() {
   // Taking points from PndEmcBump

  for (int j=0; j < fEmcInput->GetEntriesFast(); j++ ) {
    PndEmcBump* clu = (PndEmcBump*) fEmcInput->At(j);
    
    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(clu->x());  
    hit->SetY(clu->y());
    hit->SetZ(clu->z());
  
    
    if (fVerbose) cout << "EMC BUMP " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(1.);  
    hit->SetDy(1.);
    hit->SetDz(1.);

    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("EmcBump"));
    //hit->SetTrackID(point->GetTrackID());
    hit->SetTrackID(-1);

    hit->SetRefIndex(j);

    hit->SetLink(FairLink("EmcBump", j));
    
    if (fVerbose)  hit->Print();
    
  }  // end of TpcCluster loop  
} 

//_________________________________________________________________
void PndLheHitsMaker::GetGemPoints() {
   // Taking points from PndGemMCPoint
  
  if (fVerbose)
    cout << " PndLheHitsMaker::GetGemPoints(): Gem points entries " << fGemInput->GetEntriesFast() <<endl;
  
  
  for (int j=0; j < fGemInput->GetEntriesFast(); j++ ) {
    PndGemMCPoint* point = (PndGemMCPoint*) fGemInput->At(j);
    
    if (fVerbose) point->Print(" "); //PR(point->GetTrackID());

    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    if (fGemResolution<=0.)
      {
	hit->SetX(point->GetX());
	hit->SetY(point->GetY());
	hit->SetZ(point->GetZ());
	
	if (fVerbose) cout << "GEM POINT " 
			   << hit->GetX() << " " 
			   << hit->GetY() << " " 
			   << hit->GetZ() << " RADIUS " 
			   << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
	
	hit->SetDx(0.03);
	hit->SetDy(0.03);
	hit->SetDz(0.03);
      }
    else
      {
	hit->SetX(gRandom->Gaus(point->GetX(),fGemResolution));
	hit->SetY(gRandom->Gaus(point->GetY(),fGemResolution));
	hit->SetZ(gRandom->Gaus(point->GetZ(),fGemResolution));
	
	if (fVerbose) cout << "GEM POINT SMEARED " 
			   << hit->GetX() << " " 
			   << hit->GetY() << " " 
			   << hit->GetZ() << " RADIUS " 
			   << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
	
	hit->SetDx(fGemResolution);
	hit->SetDy(fGemResolution);
	hit->SetDz(fGemResolution);
      }
    
    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("GEMPoint"));
    hit->SetTrackID(point->GetTrackID());
    hit->SetRefIndex(j);

    hit->SetLink(FairLink("GEMPoint", j));
    
    if (fVerbose)  hit->Print();
  }  // end of TpcPoints loop
  
} 

//_________________________________________________________________
void PndLheHitsMaker::GetGemHits() {
   // Taking points from PndGemHits

  for (int j=0; j < fGemInput->GetEntriesFast(); j++ ) {
    PndGemHit* point = (PndGemHit*) fGemInput->At(j);
    
    PndLheHit* hit = AddHit();
    hit->SetHitNumber(fNHit++);
    
    hit->SetX(point->GetX());  
    hit->SetY(point->GetY());
    hit->SetZ(point->GetZ());
  
    
    if (fVerbose) cout << "GEM HIT " 
		       << hit->GetX() << " " 
		       << hit->GetY() << " " 
		       << hit->GetZ() << " RADIUS " 
		       << sqrt((hit->GetX()*hit->GetX())+(hit->GetY()*hit->GetY())) << "\n";
    
    hit->SetDx(point->GetDx());  
    hit->SetDy(point->GetDy());
    hit->SetDz(point->GetDz());
    
    hit->SetDetectorID(FairRootManager::Instance()->GetBranchId("GEMHit"));
    hit->SetTrackID(-1);
    if ( (fGemSimMode) && (point->GetRefIndex()!=-1) )
    {
		PndGemMCPoint* myPoint = (PndGemMCPoint*)(fGemMCArray->At(point->GetRefIndex()));
		hit->SetTrackID(myPoint->GetTrackID());
    }
    hit->SetRefIndex(j);
    hit->SetLink(FairLink("GEMHit", j));
    
    if (fVerbose)  hit->Print();
    
  }  // end of TpcCluster loop  
} 
//_________________________________________________________________
void PndLheHitsMaker::Exec(Option_t * option) {
    
  cout << "\n\n  ***  Event # " << ++fNofEvents << endl;

  cout << " =====   PndLheHitsMaker   =====\n";

  Reset();

  if ((fMvdMode==1) && (fMvdMCArray->GetEntriesFast()>0))  GetMvdPoints();
  if ((fMvdMode==2) && ((fMvdStripHitArray->GetEntriesFast()+fMvdPixelHitArray->GetEntriesFast())>0))  GetMvdHits();
  
  if ((fTpcMode==1) && (fTpcInput->GetEntriesFast()>0))  GetTpcPoints();
  if ((fTpcMode==2) && (fTpcInput->GetEntriesFast()>0))  GetTpcClusters();

  if ((fSttMode==1) && (fSttMCArray->GetEntriesFast()>0))    GetSttPoints();
  if ((fSttMode==2) && (fSttHitInput->GetEntriesFast()>0))   GetSttHit();
  if ((fSttMode==3) && (fSttHelixInput->GetEntriesFast()>0)) GetSttHelixHit();
  if ((fSttMode==4) && (fSttHelixInput->GetEntriesFast()>0)&& (fSttHitInput->GetEntriesFast()>0) ) 
    GetSttHelixHitMC();
  if ((fSttMode==5) && (fSttPRHelixInput->GetEntriesFast()>0)) GetSttPRHelixHit();
  
  if ((fEmcMode==2) && (fEmcInput->GetEntriesFast()>0))  GetEmcClusters();
  if ((fEmcMode==3) && (fEmcInput->GetEntriesFast()>0))  GetEmcBumps();
  
  if ((fGemMode==1) && (fGemInput->GetEntriesFast()>0))  GetGemPoints();
  if ((fGemMode==2) && (fGemInput->GetEntriesFast()>0))  GetGemHits();

  cout << "  Total number of hits for tracking: " <<
    setw(5) << fLheHits->GetEntriesFast() << endl;

  if (fSimulation)  CheckTracks();
  if (fVerbose) PrintTracks(0);
  
  //  fTpcInput->Clear("C");
  //  fMvdInput->Clear("C");

}

//_________________________________________________________________
PndLheHit * PndLheHitsMaker::AddHit() {
  // Creates a new hit in the TClonesArray.

  TClonesArray& hitRef = *fLheHits;
  Int_t size = hitRef.GetEntriesFast();
  return new(hitRef[size]) PndLheHit();
}

//_________________________________________________________________
void PndLheHitsMaker::Register() {
  //---

  FairRootManager::Instance()->
    Register("LheGeantTrack","Lhe", fGeantTracks, fPersistence);
  //    Register("LheGeantTrack","Lhe", fGeantTracks, kTRUE);

  FairRootManager::Instance()->
    Register("LheHit","Lhe", fLheHits, fPersistence);

}

//_________________________________________________________________
void PndLheHitsMaker::Finish() {

//   FairRootManager *fManager =FairRootManager::Instance();	
//   fManager->Fill();

  cout << " PndLheHitsMaker::Finish : Geant tracks " <<
    fGeantTracks->GetEntriesFast() << endl;

  fGeantTracks->Delete();
  fLheHits->Delete();

}

//_________________________________________________________________
void PndLheHitsMaker::Reset() {
  //---
  //  cout << " PndLheHitsMaker::Reset  " << endl;

  fNHit = 0;
  fNTrack = 0;	

  if (fLheHits->GetEntriesFast() != 0)  fLheHits->Clear("C");
  if (fGeantTracks->GetEntriesFast() != 0)fGeantTracks->Clear("C");

}

//_______________________________________________________________
void PndLheHitsMaker::CheckTracks() {

  using namespace std;

  Int_t  parent, tr_num, tr_pdg;
  Int_t nhit = 0;

  TObjArray dtracks;
  TObjArray  *chits;

  PndMCTrack *gtrack = 0;

  Int_t nMCtracks = fListMCtracks->GetEntriesFast();

  Int_t selected_tracks = 0;

  for (Int_t ih = 0; ih < fLheHits->GetEntriesFast(); ih++) {
    PndLheHit *hit = (PndLheHit *)fLheHits->UncheckedAt(ih);
    SetTrack(hit);
    //    hit->Print();
  }

  if (fVerbose) {
    cout << " MC tracks in event: " <<fListMCtracks->GetEntriesFast();
    cout << " tracks in TPC " << fNTrack << endl;
  }

  for (Int_t itrack = 0; itrack < fNTrack; itrack++) {

    PndLheCandidate *track = (PndLheCandidate*)fGeantTracks->UncheckedAt(itrack);
    tr_num = track->GetTrackNumber();

    if(tr_num >= 0 && tr_num < nMCtracks) {
      tr_pdg = 0;
      parent = 99999;

      gtrack = (PndMCTrack *) fListMCtracks->At(tr_num);

      tr_pdg = gtrack->GetPdgCode();
      track->SetPid(tr_pdg);
      track->SetCharge(TMath::Sign(1,tr_pdg));
      parent = gtrack->GetMotherID();
      if(parent == -1) {
	track->ComesFromMainVertex(kTRUE);
	track->SetGood(kTRUE);
      }

      TVector3 vertex = gtrack->GetStartVertex();
      track->SetVertex(vertex.X(), vertex.Y(), vertex.Z());

      TVector3 moment = gtrack->GetMomentum();
      track->SetPx(moment.X());
      track->SetPy(moment.Y());
      track->SetPz(moment.Z());

      nhit = track->GetNumberOfHits();
      chits = track->GetRHits();
      // print some info about track
      
      if (fVerbose>2) {
	cout << "\n *** track *** " << track->GetTrackNumber() <<
	  " pdg " << track->GetPid() <<
	  " nhits " <<nhit << " p_t ";
	printf("%f", track->GetPt());
	cout << "\n";
      }
      
      if (fVerbose) track->Print();
      PndLheHit *hit = (PndLheHit *)chits->First();
      if (fVerbose) hit->Print();

      if (fCuts->IsGoodGeantTrack(track)) {
	selected_tracks++;
      }

      if (parent == 99999 ) cout << "  track parent unknown";

    }
  }

#if 0
  for (Int_t ih = 0; ih < fLheHits->GetEntriesFast(); ih++) {
    PndLheHit *hit = (PndLheHit *)fLheHits->UncheckedAt(ih);
    Int_t mTrackNumber = hit->GetTrackID();
    for (Int_t itrack = 0; itrack < fNTrack; itrack++) {
      PndLheCandidate *track = (PndLheCandidate*)fGeantTracks->UncheckedAt(itrack);
      Int_t tr_num = track->GetTrackNumber();
      if (mTrackNumber ==  tr_num && !track->IsGood()) {
	hit->SetUsage(kTRUE);
      }
    }
  }
#endif

  cout << "Total number of tracks in TPC: " <<
    setw(5) << fGeantTracks->GetEntriesFast() << endl;
  cout << "           Good tracks in TPC: " <<
    setw(5) << selected_tracks << endl;
     
}

//________________________________________________________________
void PndLheHitsMaker::SetTrack(PndLheHit *hit) {

  Int_t mTrackNumber = hit->GetTrackID();
  Int_t mTN;

#ifdef HITS
     cout << " fNHit = " << hit->GetHitNumber();
     cout << " fNTrack = " << fNTrack;
     cout << " mTrackNumber = " << mTrackNumber;
     cout << endl;
#endif

  PndLheCandidate *track = NULL;

  Bool_t newtrack = kTRUE;
  
  if (fNTrack > 0) {
    for (Int_t i = 0; i < fGeantTracks->GetEntriesFast(); i++) {
      PndLheCandidate *currtrk = ((PndLheCandidate*) fGeantTracks->At(i));
      mTN = currtrk->GetTrackNumber();
      //      PR( currtrk->GetTrackNumber());
      if (mTN == mTrackNumber) {
	track = currtrk;
	newtrack = kFALSE;
	break;
      }
    }
  }
  if (newtrack)  track = AddTrack(mTrackNumber);

  track->AddHit(hit);
  
}

//________________________________________________________________
PndLheCandidate *PndLheHitsMaker::AddTrack(Int_t mTrackNumber) {
// Add a new track to the list of tracks for this event.

  TClonesArray &tracks = *fGeantTracks;
  PndLheCandidate *track = new(tracks[fNTrack++]) PndLheCandidate(mTrackNumber);
  return track;
}

//_________________________________________________________________
void PndLheHitsMaker::PrintTracks(Int_t ntr) {


  //  cout << " " << fNTrack << " tracks  \n\n ";
  //  cout << "Station X Y Z Px Py Pz \n\n";

  Int_t  ptracks;

  if(ntr != 0 ) 
    ptracks = ntr;
  else
    ptracks = fGeantTracks->GetEntriesFast();

  for (Int_t itrack = 0; itrack < ptracks; itrack++) {
	
    PndLheCandidate *track = (PndLheCandidate *)fGeantTracks->UncheckedAt(itrack);

    track->Print(); //

  }

}

ClassImp(PndLheHitsMaker)