#include "PndTrkPrototype.h"

#include "PndSttHit.h"
#include "PndSttPoint.h"
#include "PndSttTrack.h"
#include "PndSttHelixHit.h"
#include "PndSttTube.h"
#include "PndSttMapCreator.h"
#include "PndSttGeometryMap.h"
#include "PndTrackCand.h"
#include "PndTrackCandHit.h"
#include "PndTrack.h"
#include "PndMCTrack.h"

#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairTrackParP.h"

#include "TGeoManager.h"
#include "TClonesArray.h"
#include "TGeoVolume.h"
#include "TVector3.h"
#include "TRandom.h"
#include "TH1F.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TGeoTube.h"
#include "TF1.h"
#include "TArc.h"
#include "TPolyLine.h"
#include "TLine.h"
#include "TMarker.h"
#include "TH2F.h"
#include "TMatrixT.h"
#include "TMinuit.h"
#include "TGraph.h"

#include <iostream>
#include <sstream>
#include <string>
#include <cmath>
#include <vector>

using namespace std;

typedef vector<pair<int,int> > myvec;
typedef vector<pair<int,int> >::iterator myiter;

// ========== Default constructor ========================================
PndTrkPrototype::PndTrkPrototype() : fVerbose(0), fDisplay(kFALSE), fIsHor(kFALSE), fSmallButOk(kFALSE), fPersistence(kTRUE), FairTask("Tracking") {

  sprintf(fMCTrackBranch,  "MCTrack");
  sprintf(fSttPointBranch, "STTPoint");
  sprintf(fSttHitBranch,   "STTHit");

  // init arrays
  fSttParameters = NULL;
  fSttTubeArray  = NULL;

  // init counters
  nTotalMC = 0;
  nReco = 0;
  nRecoLess = 0;
  nRecoMore = 0;

}


// ========== Default destructor =========================================
PndTrkPrototype::~PndTrkPrototype() 
{
}


// ========== Init =======================================================
InitStatus PndTrkPrototype::Init()
{
  fEventCounter = 0;
  fNfailed = 0;

  FairRootManager * ioman = FairRootManager::Instance();
  if (!ioman) {
    cout << "-E- PndTrkPrototype::Init: "  
	 << "RootManager not instantiated, return!" << endl;
    return kFATAL;
  }

  // ---> recalling STT map
  PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters);
  fSttTubeArray = mapper->FillTubeArray();

  // get the MCTrack array
  fMCTrackArray = (TClonesArray*) ioman->GetObject(fMCTrackBranch);
  if (!fMCTrackArray)
    {
      cout << "-E- PndTrkPrototype::Init: No MCTrack array!" 
	   << endl;
      return kERROR;
    }

  // get STT points array
  fSttPointArray = (TClonesArray*) ioman->GetObject(fSttPointBranch);
  if (!fSttPointArray) {
    cout << "-E- PndTrkPrototype::Init: "
	 << "No STTPoint array!" << endl;
    return kERROR;
  }

  // get STT hits array (digi)
  fSttHitArray = (TClonesArray*) ioman->GetObject(fSttHitBranch);
  if (!fSttHitArray) {
    cout << "-E- PndTrkPrototype::Init: "
	 << "No STTHit array, return!" << endl;
    return kERROR;
  }

  // definition of output histograms
  hdist        = new TH1F("hdist",       "hdist",       100, 0, 5);
  hresfit      = new TH1F("hresfit",     "hresfit",     200, -3., 3.);
  hresrefit    = new TH1F("hresrefit",   "hresrefit",   200, -3., 3.);
  hlastresid   = new TH1F("hlastresid",  "hlastresid",  200, -3., 3.);
  chi2r_prefit = new TH1F("chi2r_prefit","chi2r_prefit",200, 0, 20);
  chi2r_fit    = new TH1F("chi2r_fit",   "chi2r_fit",   200, 0, 20);
  chi2r_refit  = new TH1F("chi2r_refit", "chi2r_refit", 200, 0, 20);
  chi2r_oldfit = new TH1F("chi2r_oldfit","chi2r_oldfit",200, 0, 20);
  cumul_fit    = new TH1F("cumul_fit",   "cumul_fit",   200, 0, 1);
  cumul_refit  = new TH1F("cumul_refit", "cumul_refit", 200, 0, 1);
  cumul_oldfit = new TH1F("cumul_oldfit","cumul_oldfit",200, 0, 1);
  hstats       = new TH1F("hstats",      "hstats",      5,  -0.5, 4.5);
  hmeanres     = new TH1F("hmeanres",    "hmeanres",    200, 0, 1);
  htrklength   = new TH1F("htrklength",  "htrklength",  200, 0, 18);
  hde          = new TH1F("hde",         "hde",         200, 0., 500);
  hdedx        = new TH1F("hdedx",       "hdedx",       200, 0., 50);
  htubemult    = new TH1F("htubemult",   "htubemult",   128, 0.5, 128.5);
  hgoodd       = new TH1F("hgoodd",      "hgoodd",      100, 0, 0.5);

  // eventually... create and register output (TrackCand or Track)
  // .... it goes here!

  return kSUCCESS;

}

// ========== Retrieve params ============================================
void PndTrkPrototype::SetParContainers() {
  FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb();
  fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
}


// ========== Reconstruct ================================================
void PndTrkPrototype::Exec(Option_t* opt) {

  nPrimaryMCtrks = 0;
  nRecoTrks = 0;
  fAdjacent = kFALSE;
  fSmallButOk = kFALSE;

  for (int i=0; i<fNTubes+1; i++) {
    for (int ii=0; ii<fNTubes+1; ii++) {
      IsAdded[i][ii]   = kFALSE;
      IsChecked[i][ii] = kFALSE;
      IsFake[i][ii]    = kFALSE;
    }
    IsUsed[i] = kFALSE;
  }
  for (int i=0; i<fNMAXCLUSTERS+1; i++)
    fIsPrimary[i] = kFALSE;

  if ((fVerbose==0 && fEventCounter%1000==0) || (fVerbose>0))
    cout << "=============== Event # " << fEventCounter << endl;

  if (fDisplay) {
    if (fEventCounter==0)
      c = new TCanvas("display_zy","display_zy", 0, 0, 1000, 1000);
    DrawGeometry();
  }

  fEventCounter++;

  PndMCTrack* pMCtr;
  PndSttPoint*pSttPnt;
  PndSttHit *pSttHit;
  PndSttTube *pSttTube;
  myvec hitlist;
  myvec middlelist;
  myvec lowerlist;
  myvec upperlist;

  Double_t time;

  Int_t nSttMCPoint = fSttPointArray->GetEntriesFast();
  Int_t nMCtracks = fMCTrackArray->GetEntriesFast();
  if (fVerbose > 0 && nMCtracks > 0) {
    cout << "# of MC tracks " << nMCtracks << endl;
    cout << "# of STT MC points: " << nSttMCPoint << endl;
  }

  Double_t MCpoints[nMCtracks][nSttMCPoint][2];

  // fetching the STT hits
  fNSttHits = fSttHitArray->GetEntriesFast();
  if (fVerbose>0)
    cout << "# of STT hits: " << fNSttHits << endl;
  if (fNSttHits > fNMAXSIZE) {
    if (fVerbose>1)
      cout << "There are too many hits (" << fNSttHits << ") in this event... --> discarded!" << endl;
    return;
  }
  if (fNSttHits==0)
    return;

  Int_t nMCpt[nMCtracks];
  for (Int_t imc=0; imc<nMCtracks; imc++) 
    nMCpt[imc] = 0;

  for (Int_t ipt=0; ipt<nSttMCPoint; ipt++) {
    pSttPnt = (PndSttPoint*) fSttPointArray->At(ipt);
    Int_t trackID = pSttPnt->GetTrackID();
    TVector3 PntCoord;
    pSttPnt->Position(PntCoord);
    MCpoints[trackID][nMCpt[trackID]][0] = PntCoord.Z();
    MCpoints[trackID][nMCpt[trackID]][1] = PntCoord.Y();
    nMCpt[trackID]++;
  }

  Double_t positions[2][fNMAXSIZE];
  TMatrixT<Double_t> hitArray(fNSttHits,8);
  TMatrixT<Double_t> MiddleArray;
  for (Int_t imc=0; imc<nMCtracks; imc++) {
    pMCtr = (PndMCTrack*) fMCTrackArray->At(imc);
    //    if (pMCtr->GetMotherID()==-1 &&  nMCpt[imc]!=0)
    if (pMCtr->GetMotherID()==-1 && nMCpt[imc]>2)
      nPrimaryMCtrks++;
  }

  Int_t ylayer[fNLayers+1][2];
  for (int iy=0; iy<fNLayers+1; iy++) {
    ylayer[iy][0] = 0; ylayer[iy][1] = 0; }
  Int_t HowMany = 1;
  Bool_t yFound = kFALSE;
  Int_t layerID = 0;
  vector<int> myvector;
  myvec fakehits;

  for (Int_t i=0; i<fNSttHits; i++) {
    
    pSttHit = (PndSttHit*) fSttHitArray->At(i);
    if (!pSttHit) continue;

    yFound = kFALSE;

    // hit info
    Double_t isochrone = pSttHit->GetIsochrone();
    Int_t tubeID       = pSttHit->GetTubeID();
    hitlist.push_back(make_pair(tubeID,i));
    htubemult->Fill(tubeID);

    pSttTube = (PndSttTube*) fSttTubeArray->At(tubeID);
    layerID = pSttTube->GetLayerID();
    TVector3 center_pos = pSttTube->GetPosition();
    positions[0][i] = center_pos.Z();
    positions[1][i] = center_pos.Y();

    hitArray[i][0] = tubeID;
    hitArray[i][1] = i;
    hitArray[i][2] = center_pos.Z();
    hitArray[i][3] = center_pos.Y();
    hitArray[i][4] = pSttHit->GetIsochrone();
    hitArray[i][5] = pSttHit->GetIsochroneError();
    hitArray[i][6] = pSttHit->GetEnergyLoss();
    hitArray[i][7] = pSttHit->GetRefIndex();

    if (i==0) {
      ylayer[i][0] = layerID;
      ylayer[i][1]++;
      myvector.push_back(layerID);
    }
    else {
      for (int iy=0; iy<HowMany; iy++) {
	if (layerID == ylayer[iy][0]) {
	  ylayer[iy][1]++;
	  yFound = kTRUE;
	  break;
	}
      }
      if (yFound==kFALSE) {
	ylayer[HowMany][0] = layerID;
	ylayer[HowMany][1]++;
	HowMany++;
	myvector.push_back(layerID);
      }
    }
  
    if (fVerbose>3) 
      cout << "Hit # " << i << "\t tube " << tubeID << "\t pos: (" << center_pos.Z() 
 	   << ", " << center_pos.Y() << ")" << "\t isochrone " << isochrone << endl;
  }

  if (fDisplay)
    DrawIsochrone(hitArray, fNSttHits, -1, 0);

  std::sort(myvector.begin(), myvector.end());

  if (fVerbose>1) {
    cout << endl;
    cout << "==== Map before hit ordering [tubeID, hit#]" << endl;
    for (myiter it = hitlist.begin(); it!= hitlist.end(); it++)
      cout << "(" << it->first << ", " << it->second << ")  ";
    cout << endl; cout << endl;
    if (HowMany==1)
      cout << "1 layer has been hit" << endl;
    else 
      cout << HowMany << " layers have been hit" << endl;
    for (int iy=0; iy<HowMany; iy++) 
      cout << "Layer with y = " << ylayer[iy][0] << " has been hit " << ylayer[iy][1] << " times" << endl;
    cout << endl;
  }

  vector<int>::iterator itl; 
  vector<int>::iterator itm;
  vector<int>::iterator itu;
  if (HowMany==3 && fNSttHits>30) {
    itl = myvector.begin();
    itm = myvector.begin()+1;
    itu = myvector.end()-1;

    if ((*itm-*itl==1) && (*itu-*itm==1)) {
      fAdjacent = kTRUE;
      for (int i=0; i<fNSttHits; i++) {
	pSttTube = (PndSttTube*) fSttTubeArray->At(hitArray[i][0]);
	layerID = pSttTube->GetLayerID();
	if (layerID==*itl || layerID==*itu) {
	  fakehits.push_back(make_pair(hitArray[i][0], hitArray[i][1]));
	  IsFake[(int)hitArray[i][0]][(int)hitArray[i][1]] = kTRUE;
	}
	else if (layerID==*itm) {
	  for (myiter ii=middlelist.begin(); ii!=middlelist.end(); ii++) {
	    if (hitArray[i][0] == ii->first) {
	      fakehits.push_back(make_pair(hitArray[i][0], hitArray[i][1]));
	      IsFake[(int)hitArray[i][0]][(int)hitArray[i][1]] = kTRUE;
	      break;
	    }
	  }
	  if (IsFake[(int)hitArray[i][0]][(int)hitArray[i][1]] == kFALSE) {
	    middlelist.push_back(make_pair(hitArray[i][0], hitArray[i][1]));
	    Int_t arraypos = middlelist.size();
	    MiddleArray.ResizeTo(arraypos,8);
	    for (int k=0; k<8; k++) 
	      MiddleArray[arraypos-1][k] = hitArray[i][k];
	  }
	}
      }
    }
  }

  TArrayI list_multi;
  Int_t size_multi = 0;
  Bool_t IsInMultiList = kFALSE;
  for (Int_t i=0; i<fNSttHits; i++) {
    for (Int_t j=i+1; j<fNSttHits; j++) {
      if (hitArray[i][0] == hitArray[j][0]) {
	if (fAdjacent==kTRUE) {
	  fakehits.push_back(make_pair(hitArray[i][0], hitArray[i][1]));
	  break;
	}
 	size_multi = list_multi.GetSize();
 	list_multi.Set(size_multi+1);
	IsInMultiList = kFALSE;
	for (Int_t ii=0; ii<list_multi.GetSize(); ii++) {
	  if (hitArray[i][0] == list_multi.At(ii)) {
	    IsInMultiList = kTRUE;
	    break;
	  }
	}
	if (IsInMultiList == kFALSE) {
	  size_multi = list_multi.GetSize();
	  list_multi.Set(size_multi+1);
	  list_multi.AddAt(hitArray[i][0],size_multi);
	}
      }
    }
  }
  if (fVerbose>1) {
    if (list_multi.GetSize()!=0) {
      cout << "List of tubes with multi hits: " << endl;
      for (Int_t ii=0; ii<list_multi.GetSize(); ii++)
	cout << list_multi.At(ii) << "\t" ;
      cout << endl;
    }
  }

  // MC info used only to display MC tracks (for comparison with the results)
  Double_t mMC[nMCtracks], qMC[nMCtracks];
  Bool_t fMC[nMCtracks];
  TMatrixT<Double_t> MCpts;
  for (Int_t imc=0; imc<nMCtracks; imc++) {
    fMC[imc] = kFALSE;
    mMC[imc] = 0.;   
    qMC[imc] = 0.;
    if (nMCpt[imc]>0) {
      MCpts.ResizeTo(nMCpt[imc], 2);
      for (Int_t ik=0; ik<nMCpt[imc]; ik++) {
	MCpts[ik][0] = MCpoints[imc][ik][0];
	MCpts[ik][1] = MCpoints[imc][ik][1];
	if (fVerbose>3) cout << "MC point " << ik << ": (" << MCpts[ik][0] << "," << MCpts[ik][1] << ")" << endl;
      }
      if ((MCpoints[imc][1][0] - MCpoints[imc][0][0]) != 0) {
	fMC[imc] = kTRUE;
	mMC[imc] = (MCpoints[imc][nMCpt[imc]-1][1] - MCpoints[imc][0][1]) / (MCpoints[imc][nMCpt[imc]-1][0] - MCpoints[imc][0][0]);
	qMC[imc] = 10*MCpoints[imc][0][1] - mMC[imc]*10*MCpoints[imc][0][0];
      }
      if (fDisplay && fMC[imc]) 
	DrawMCPoint(MCpts,nMCpt[imc], mMC[imc], qMC[imc]);
    }
  }

  myvec orderedMap[fNMAXCLUSTERS];
  Int_t sizecl[fNMAXCLUSTERS];
  for (Int_t i=0; i<fNMAXCLUSTERS; i++)     
    sizecl[i] = 0;
  Int_t clidx = 0;
  if (fAdjacent==kFALSE)
    clidx = CreatingClusters(fNSttHits, hitArray, list_multi, orderedMap, sizecl, fakehits);
  else 
    clidx = CreatingClusters(middlelist.size(), MiddleArray, list_multi, orderedMap, sizecl, fakehits);
  if (fVerbose>1) {
    cout << endl;
    cout << "==== Map after hit ordering and clustering (tubeID, hit#)" << endl;
    if (clidx==1)
      cout << "There is " << clidx << " cluster" << endl;
    else
      cout << "There are " << clidx << " clusters" << endl;
    for (Int_t k=0; k<clidx; k++) {
      cout << "-- Cluster # " << k << " --> size = " << orderedMap[k].size() << endl;
      for (myiter itt=orderedMap[k].begin(); itt!=orderedMap[k].end(); itt++) 
	cout << "(" << itt->first << ", " << itt->second << ")  " ;
    }
    cout << endl;     cout << endl;
  }

  Int_t size = 0;
  for (Int_t icl=0; icl<clidx; icl++) {
    if (fVerbose>1) {
      cout << "*** Total number of clusters: " << clidx << endl; cout << endl; }
    if (icl>=10) {
      if (fVerbose>2)
	cout << "There are too many (" << clidx << ") clusters ... discard event" << endl;
      break;
    }

    if (fVerbose>1)
      cout << "*** Reconstructing cluster " << icl << " ***" << endl;

    fSmallButOk = kFALSE;

    size = sizecl[icl];

    // if the cluster is small, check if the hits are at the border
    // if so, the tracking is meaningful
    if (size>2 && size<5) {
      Int_t checkBorder = 0;
      Int_t nBorder = 0;
      for (myiter it = orderedMap[icl].begin(); it!=orderedMap[icl].end(); it++) {
	pSttTube = (PndSttTube*) fSttTubeArray->At(it->first);
	checkBorder = pSttTube->IsLayerLimit();
	if (checkBorder==1) nBorder++;
      }
      if (nBorder==size)
	fSmallButOk = kTRUE;
    }

    myvec cluster;
    Int_t clsize = 0;
    TMatrixT<Double_t> hit; 
    cluster = orderedMap[icl];
    clsize = cluster.size();
    if (fVerbose>1)
      cout << "Cluster " << icl << " \t size: " << clsize << endl;
    // Filling array with hit positions
    hit.ResizeTo(clsize,8);
    Int_t j = 0;
    for (myiter it=cluster.begin(); it!=cluster.end(); it++) {
      hit[j][0] = it->first;  // tubeID
      hit[j][1] = it->second; // hit #
      hit[j][2] = ( (PndSttTube*) fSttTubeArray->At(it->first))->GetPosition().Z();
      hit[j][3] = ( (PndSttTube*) fSttTubeArray->At(it->first))->GetPosition().Y();
      hit[j][4] = ( (PndSttHit*) fSttHitArray->At(it->second) )->GetIsochrone();
      hit[j][5] = ( (PndSttHit*) fSttHitArray->At(it->second) )->GetIsochroneError();
      hit[j][6] = ( (PndSttHit*) fSttHitArray->At(it->second) )->GetEnergyLoss();
      hit[j][7] = ( (PndSttHit*) fSttHitArray->At(it->second) )->GetRefIndex();
  
      if (fVerbose>3) 
	cout << "Coordinates of hit " << hit[j][1] << " [tube " << hit[j][0] << "]: (" 
	     << hit[j][2] << ", " << hit[j][3] << ")\t isochrone = (" << hit[j][4] 
	     << " +/- " << hit[j][5] << ")" << endl;
      j++;
    }
    
    if (fDisplay) 
      DrawIsochrone(hit, clsize, 0, icl);

    // ------- PREFIT using the centers of the firing tubes
    Double_t ppref[2], pminfit[2], pminfit2[2], pminfit3[2];
    Double_t errpminfit[2], errpminfit2[2], errpminfit3[2];
    Double_t pfinal[2], errpfinal[2];
    for (int i=0; i<2; i++) {
      ppref[i] = 0.; pminfit[i] = 0.; 
      pminfit2[i] = 0.; pminfit3[i] = 0.; 
      pfinal[i] = 0.; errpminfit2[i] = 0.;
      errpminfit3[i] = 0.; errpfinal[i] = 0.; 
    }

    Double_t mdp;
    Double_t residual = 0., sumres = 0., resmed = 0.;
    Double_t chisq = 0.;
    Bool_t fPrefit, fCheckD, fCheckR, fMinFit, fIntFind, fRefit;
    Int_t iter = 0;
    Int_t fake = 0;
    Int_t fcnversion = 0;
    TMatrixT<Double_t> hitnew(clsize,8);

    fPrefit  = kFALSE;
    fCheckD  = kFALSE;
    fCheckR  = kFALSE;
    fMinFit  = kFALSE;
    fIntFind = kFALSE;
    fRefit   = kFALSE;

    do {

      fPrefit = LinFitPerp(hit, clsize, ppref, mdp);

      if (fPrefit) {
	// 	if (iter == 0 && fDisplay) 
	// 	  DrawFit(1, iter, ppref, icl);
	fCheckD = CheckDistance(clsize, hit, ppref, fake, hitnew, fakehits);
	if (fCheckD) {
	  iter = 1;
	  clsize -= fake;
	  hit.ResizeTo(clsize,8);
	  hit = hitnew;
	  if (fVerbose>2) {
	    cout << "Size after checkdistance " << clsize << endl;
	    for (Int_t gg=0; gg<clsize; gg++)
	      cout << hitnew[gg][0] << "\t";
	    cout << endl;
	  }
	}
	else 
	  fCheckD = kFALSE;
      }
      else 
	fCheckD = kFALSE;
    } while (fCheckD);

    //     if (fDisplay) {
    //       DrawIsochrone(hit, clsize, iter, icl);
    //     if (fPrefit)
    //       DrawFit(1, iter, ppref, icl);
    //    }

    iter = 0;
    if (fPrefit) {
      do {
	fcnversion = 1;
	fMinFit = MinuitFit(fcnversion, clsize, hit, ppref, pminfit, errpminfit, chisq);
	if (fMinFit) {
	  //	  if (iter==0 && fDisplay) 
	  //	   	  DrawFit(2, iter, pminfit, icl);
	  fCheckR = CheckResiduals(2, clsize, hit, pminfit, fake, hitnew, fakehits);
	  if (fCheckR) {
	    iter = 1;
	    clsize -= fake;
	    hit.ResizeTo(clsize,8);
	    if (clsize>0)
	      hit = hitnew;
	    if (fVerbose>2) { 
	      cout << "Size after checkresiduals " << clsize << endl;
	      for (Int_t gg=0; gg<clsize; gg++)
		cout << hitnew[gg][0] << "\t";
	      cout << endl;
	    }
	  }
	}
	else 
	  fCheckR = kFALSE;
      } while (fCheckR);

      iter = 0;
      TMatrixT<Double_t> inters(clsize,8);
      fIntFind = IntersectionFinder(clsize, hit, pminfit, inters);
      if (fIntFind) {
	// 	if (fDisplay)
	// 	  DrawIntersections(inters, clsize, iter, icl);
	fcnversion = 2;
	fRefit = MinuitFit(fcnversion, clsize, inters, pminfit, pminfit2, errpminfit2, chisq);
	if (!fRefit) {
	  pminfit2[0] = pminfit[0]; pminfit2[1] = pminfit[1]; 
	  errpminfit2[0] = errpminfit[0]; errpminfit2[1] = errpminfit[1];				  
	  fRefit = kTRUE;
	}
	if (fRefit) {
	  // 	if (fDisplay)
	  // 	  DrawFit(3, iter, pminfit2, icl);
	  if (chisq>5.) {
	    fcnversion = 3;
	    fRefit = MinuitFit(fcnversion, clsize, inters, pminfit, pminfit2, errpminfit2, chisq);
	    
	    if (!fRefit) {
	      pminfit2[0] = pminfit[0]; pminfit2[1] = pminfit[1]; 
	      errpminfit2[0] = errpminfit[0]; errpminfit2[1] = errpminfit[1];
	      
	      fRefit = kTRUE;
	    }
	  }

	  iter = 1;
	  fIntFind = IntersectionFinder(clsize, hit, pminfit2, inters);
	  if (fIntFind) {
	    // 	  if (fDisplay)
	    // 	    DrawIntersections(inters, clsize, iter, icl);
	    fcnversion = 4;
	    fRefit = MinuitFit(fcnversion, clsize, inters, pminfit2, pminfit3, errpminfit3, chisq);
	    if (!fRefit) {
	      pminfit3[0] = pminfit2[0]; pminfit3[1] = pminfit2[1];
	      errpminfit3[0] = errpminfit2[0]; errpminfit3[1] = errpminfit2[1];
	      fRefit = kTRUE;
	    }

	    if (fRefit && clsize>0) {
	      Double_t tracklength = 0., tracklengthsymm = 0.;
	      Double_t eloss = 0., elosssymm = 0.;
 	      Double_t dedx = 0.;
	      TMatrixT<Double_t> newhit, newhitsymm;
	      Int_t newsize = 0, newsizesymm = 0;
	      Double_t pRefit[2], errpRefit[2];
	      Double_t pRefSymm[2], errpRefSymm[2];
	      Double_t chisqrefit = 999.;
	      Double_t chisqrefsymm = 999.;

	      Int_t nextra = TrackIntersections(fNSttHits, pminfit3, hitArray, hit, newhit, clsize, newsize, tracklength, eloss);
	      if (fVerbose>2)
		cout << "New size after TrackIntersections(): " << newsize << " (it was " << clsize << ") " << endl;
	      if (clsize<newsize) {

		pRefit[0] = 0.; pRefit[1] = 0.;
		errpRefit[0] = 0.; errpRefit[1] = 0.;

		if (fVerbose>2) {
		  cout << "%%%%%%%%%%%%%%%%%%%%%%%" << endl;
		  cout << "Fitting new array with " << newsize << " tubes:" << endl;
		  for (int jj=0; jj<newsize; jj++)
		    cout << newhit[jj][0] << "\t";
		  cout << endl;
		}

		inters.ResizeTo(newsize, 8);
		fIntFind = IntersectionFinder(newsize, newhit, pminfit3, inters);
		Bool_t refit = MinuitFit(fcnversion, newsize, inters, pminfit3, pRefit, errpRefit, chisqrefit);
		if (refit) {
		  if (fVerbose>2) {
		    cout << "INPUT PARAMETERS" << endl;
		    cout << pminfit3[0] << ", " << pminfit3[1] << ", chi2 " << chisq << endl;
		    cout << "OUTPUT PARAMETERS (with new array) " << endl;
		    cout << pRefit[0] << ", " << pRefit[1] << ", chi2 refit " << chisqrefit << endl;
		    cout << "chisq diff " << chisqrefit - chisq << endl;
		    cout << "%%%%%%%%%%%%%%%%%%%%%%" << endl;
		  }
		}
	      }

	      Double_t psymm[2];
	      psymm[0] =  2*hit[0][3] - pminfit3[0];
	      psymm[1] = -pminfit3[1];
	      if (fVerbose>2)
		cout << "Checking intersections for the symmetric track (wrt the layer axis): " << psymm[0] << ", " << psymm[1] << endl;
	      Int_t nextrasymm = TrackIntersections(fNSttHits, psymm, hitArray, hit, newhitsymm, clsize, newsizesymm, tracklengthsymm, elosssymm);
	      if (fVerbose>2)
		cout << "New size after TrackIntersections() with symmetric parameters: " << newsizesymm << endl;
	      if (clsize<newsizesymm) {

		pRefSymm[0] = 0.; pRefSymm[1] = 0.;
		errpRefSymm[0] = 0.; errpRefSymm[1] = 0.;

		if (fVerbose>1) {
		  cout << "%%%%%%%%%%%%%%%%%%%%%%%" << endl;
		  cout << "Fitting new array with " << newsizesymm << " tubes:" << endl;
		  for (int jj=0; jj<newsizesymm; jj++)
		    cout << newhitsymm[jj][0] << "\t";
		  cout << endl;
		}
		inters.ResizeTo(newsizesymm, 8);
		fIntFind = IntersectionFinder(newsizesymm, newhitsymm, psymm, inters);
		Bool_t refitsymm = MinuitFit(fcnversion, newsizesymm, inters, psymm, pRefSymm, errpRefSymm, chisqrefsymm);
		if (refitsymm) {
		  if (fVerbose>2) {
		    cout << "INPUT PARAMETERS" << endl;
		    cout << psymm[0] << ", " << psymm[1] << ", chi2 " << chisq << endl;
		    cout << "OUTPUT PARAMETERS (with new arrya) " << endl;
		    cout << pRefSymm[0] << ", " << pRefSymm[1] << ", chi2 refit " << chisqrefsymm << endl;
		    cout << "chisq diff " << chisqrefsymm - chisq << endl;
		    cout << "%%%%%%%%%%%%%%%%%%%%%%" << endl;
		  }
		}
	      }

	      Bool_t pflag = 0;
	      if (chisqrefit<chisqrefsymm) {
		pfinal[0] = pRefit[0]; pfinal[1] = pRefit[1];
		errpfinal[0] = errpRefit[0]; errpfinal[1] - errpRefit[1];
		pflag = 1;
		hit.ResizeTo(newsize, 8);
		hit = newhit;
		clsize = newsize;
	      }
	      else if (chisqrefsymm<chisqrefit) {
		pfinal[0] = pRefSymm[0]; pfinal[1] = pRefSymm[1];
		errpfinal[0] = errpRefSymm[0]; errpfinal[1] = errpRefSymm[1];
		pflag = 1;
		hit.ResizeTo(newsizesymm, 8);
		hit = newhitsymm;
		clsize = newsizesymm;
	      }
	      else {
		if (nextra > nextrasymm) {
		  pfinal[0] = psymm[0]; pfinal[1] = psymm[1];
		}
		else {
		  pfinal[0] = pminfit3[0]; pfinal[1] = pminfit3[1];
		  errpfinal[0] = errpminfit3[0]; errpfinal[1] = errpminfit3[1];
		}
	      }
	      Int_t nfinal = TrackIntersections(fNSttHits, pfinal, hitArray, hit, newhit, clsize, newsize, tracklength, eloss);	      for (int ns=0; ns<newsize; ns++) {
		IsAdded[(int)newhit[ns][0]][(int)newhit[ns][1]] = kTRUE;
		IsUsed[(int)newhit[ns][1]] = kTRUE;
	      }

	      if (fVerbose>2) {
		cout << "Final cluster (after last TrackIntersections()): " << nfinal << endl;
		for (int gg=0; gg<newsize; gg++) 
		  cout << newhit[gg][0] << "\t";
		cout << endl;
	      }

	      if (fDisplay)
		DrawIsochrone(newhit, newsize, 1, icl);

	      // dE/dx calculation
	      htrklength->Fill(tracklength);
	      hde->Fill(eloss);
	      if (tracklength!=0) {
		dedx = eloss/tracklength;
		hdedx->Fill(dedx);
	      }

	      // filing residual histo
	      fake = 0;  
	      fCheckR = CheckResiduals(3, clsize, hit, pfinal, fake, hitnew, fakehits);
	      for (Int_t kk=0; kk<clsize; kk++) {
		residual = Point2Line(hit[kk][2], hit[kk][3], pfinal) - hit[kk][4];
		hlastresid->Fill(residual*10);
		if (residual<0.1) {
		  hgoodd->Fill(Point2Line(hit[kk][2], hit[kk][3], pfinal));
		  //		  time = GetTimeInfo(hit[kk][4]);
		  // 		  cout << "Radius " << hit[kk][4] << "\tTime " << GetTimeInfo(hit[kk][4]) << endl;
		}
		if (fVerbose>3) cout << "residual " << residual << endl;
		sumres += TMath::Abs(residual);
	      }
	      resmed = sumres/clsize;
	      hmeanres->Fill(resmed*10);//mm
	    }
	  }
	}
      }
    }
    if (fVerbose>1)
      cout << "*** Final track parameters (y=ax+b): a = " << pfinal[1] << ", b = " << pfinal[0] << " *** " << endl;
    if (fDisplay && fRefit)
      DrawFit(3, 2, pfinal, icl);

    if (fRefit)
      nRecoTrks++;

    Bool_t reverse;
    if (fVerbose>0) {
      cout << "*************** RESULTS for track " << icl << " ***************" << endl;
      if (!fMC)
	cout << "Sorry, no hits in this event... " << endl;
      else {
	if (nMCtracks>0)
	  cout << "MC values  ----> a = " << mMC[icl] << ", b = " << qMC[icl] << endl;
	if (fRefit) {
	  reverse = kFALSE;
	  cout << "RECO values ---> a = " << pfinal[1] << ", b = " << pfinal[0] << endl;
	  if (((mMC[icl]>0 && pfinal[1]<0) || (mMC[icl]<0 && pfinal[1]>0)) && clsize>0) {
	    reverse = kTRUE;
	    cout << "RECO values --> a = " << -pfinal[1] << ", b = " << 2*hit[0][3] - pfinal[0] << endl;
	  }
	}
	else
	  cout << "The track has not been reconstructed " << endl;
	cout << "***************************************" << endl;
      }
      cout << endl;
    }

    if (fVerbose>2) {
      cout << "There were " << fakehits.size() << " fake hits out of " << sizecl[icl] << " --- tubes (hit#): " << endl;
      for (myiter iitt=fakehits.begin(); iitt!=fakehits.end(); iitt++)
	cout << iitt->first << "(" << iitt->second << ")\t";
      cout << endl;
    }
    if (fakehits.size()>0 && fakehits.size()!= sizecl[icl]) {
      TMatrixT<Double_t> fakeArray(fakehits.size(),8);
      Int_t counter = 0;
      if (fVerbose>2)	{
	cout << "Creating array of fake hits " << endl;
	cout << "Tubes (hit#): " << endl;
	for (myiter ifake=fakehits.begin(); ifake!=fakehits.end(); ifake++) 
	  cout << ifake->first << "(" << ifake->second << ")\t";
	cout << endl;
      }
      for (myiter ifake=fakehits.begin(); ifake!=fakehits.end(); ifake++) {
	pSttTube = (PndSttTube*) fSttTubeArray->At(ifake->first);
	TVector3 center_pos = pSttTube->GetPosition();
	fakeArray[counter][0] = ifake->first;
	fakeArray[counter][1] = ifake->second;
	fakeArray[counter][2] = center_pos.Z();
	fakeArray[counter][3] = center_pos.Y();
	fakeArray[counter][4] = ( (PndSttHit*) fSttHitArray->At(ifake->second) )->GetIsochrone();
	fakeArray[counter][5] = ( (PndSttHit*) fSttHitArray->At(ifake->second) )->GetIsochroneError();
	fakeArray[counter][6] = ( (PndSttHit*) fSttHitArray->At(ifake->second) )->GetEnergyLoss();
	fakeArray[counter][7] = ( (PndSttHit*) fSttHitArray->At(ifake->second) )->GetRefIndex();
	counter++;
      }

      myvec newOrdMap[fNMAXCLUSTERS];
      Int_t newSize[fNMAXCLUSTERS];
      Int_t newcl = CreatingClusters(fakehits.size(), fakeArray, list_multi, newOrdMap, newSize, fakehits);
      for (Int_t nfakecl=0; nfakecl<newcl; nfakecl++) {
	orderedMap[clidx+nfakecl] = newOrdMap[nfakecl];
	sizecl[clidx+nfakecl] = newSize[nfakecl];
      }
      if (fVerbose>1) 
	if (newcl==1)
	  cout << newcl << " new cluster has been added" << endl;
	else if (newcl>1)
	  cout << newcl << " new clusters have been added" << endl;
      clidx += newcl;
    }
  }


//   hitArray.Clear();
//   MiddleArray.Clear();


  if (fVerbose>0) {
    cout << "=============== SUMMARY for event " << fEventCounter-1 << ":" << endl;
    cout << "There are " << nPrimaryMCtrks << " primary MC tracks and " << nRecoTrks << " reco tracks " << endl;
    cout << "======================================================================" << endl;
    cout << endl;
  }

  if (nPrimaryMCtrks>0) {
    nTotalMC++;
    if (nPrimaryMCtrks==nRecoTrks)
      nReco++;
    else if (nRecoTrks<nPrimaryMCtrks) {
      nRecoLess++;
      if (fVerbose>1)
	cout << "@@@@ LESS!!!! - Event " << fEventCounter-1 << "\tMC " << nPrimaryMCtrks << ", reco " << nRecoTrks << ", nSttPoints " << nSttMCPoint << endl;
    }
    else if (nRecoTrks>nPrimaryMCtrks) {
      nRecoMore++;
      if (fVerbose>1)
	cout << "#### MORE!!!! - Event " << fEventCounter-1 << "\tMC primary " << nPrimaryMCtrks 
	   << ", MC total " << nMCtracks << ", reco " << nRecoTrks << endl;
    }
  }

  // at the end of exec
  if (fDisplay) {
    char goOnChar;
    cout << "Please, press any key to show next event [type 'q' to exit]: " << endl;
    cin >> goOnChar;
    h->Delete();
    if (goOnChar == 'q') 
      exit(0);
  }

} // end exec


// ========== Method to create hit cluster ===============================
Int_t PndTrkPrototype::CreatingClusters(Int_t endloop, TMatrixT<Double_t> hit, TArrayI listmulti, vector<pair<int,int> > ordMap[fNMAXCLUSTERS+1], Int_t size[fNMAXCLUSTERS+1], vector<pair<int,int> > &fakelist) {

  if (fVerbose>1) 
    cout << "=== PndTrkPrototype::CreatingClusters " << endl;

  Int_t ncl = 0;
  Double_t min = 1000., miny = -1000.;
  Int_t ind = 0;
  Int_t itube = 0;
  PndSttTube* current_tube;
  Int_t hits_left = endloop;

  for (int i=0; i<fNMAXCLUSTERS+1; i++)
    size[i] = 0;

  do {
    fIsPrimary[ncl] = kFALSE;
    min = 1000;
    miny = -1000;
    ind = fNMAXSIZE+1;
    itube = 0; // there is no tube numbered 0

    for (int i=0; i<endloop; i++) {
      if (fVerbose>2)
 	cout  << "Hit " << hit[i][1] << "\t tube " << hit[i][0] << "\t used? " << IsAdded[(int)hit[i][0]][(int)hit[i][1]] << endl;
      if (IsAdded[(int)hit[i][0]][(int)hit[i][1]] == kTRUE)
	continue;
      current_tube = (PndSttTube*) fSttTubeArray->At(hit[i][0]);
      if (current_tube->IsSectorLimit() == -1) {
	if (hit[i][2]<min) {
	  ind = hit[i][1];
	  itube = hit[i][0];
	  min = hit[i][2];
	  miny = hit[i][3];
	}
	else if (hit[i][2]==min && hit[i][3]<miny) {
	  ind = hit[i][1];
	  itube = hit[i][0];
	  min = hit[i][2];
	  miny = hit[i][3];
	}
      }
    }

    if (itube!=0 && ind!=fNMAXSIZE)
      fIsPrimary[ncl] = kTRUE;
    else 
      fIsPrimary[ncl] = kFALSE;
    if (fVerbose>2) cout << "Is the track primary? " << fIsPrimary[ncl] << endl;
    if (fIsPrimary[ncl]==kFALSE) {
      hits_left = 0;
      ncl--;
      continue;
    }

    ordMap[ncl].push_back(make_pair(itube, ind));
    if (fVerbose>2) 
      cout << "Cluster begins with tube " << itube << " (hit " << ind << ")" << endl;
    IsAdded[itube][ind] = kTRUE;
    IsFake[itube][ind] = kFALSE;
    Bool_t Found = kFALSE;
    Int_t index;

    for (myiter jt=ordMap[ncl].begin(); jt!=ordMap[ncl].end(); jt++) {
      if (fVerbose>3)
	cout << "Tube " << jt->first << " is checked? " << IsChecked[jt->first][jt->second] << endl;
      if (IsChecked[jt->first][jt->second] == kTRUE) continue;
      itube = jt->first;
      index = jt->second;
      if (fVerbose>3)
	cout << "**** Going through the neighbours of the last added tube: tube " << jt->first << endl;
      current_tube = (PndSttTube*) fSttTubeArray->At(jt->first);
      TArrayI neighborings = current_tube->GetNeighborings();
      for (int ineigh=0; ineigh<neighborings.GetSize(); ineigh++) {
	if (fVerbose>3) 
	  cout << "---> Checking tube " << neighborings.At(ineigh) << ": has it been hit? " << endl;
	for (int ii=0; ii<endloop; ii++) {
	  if (fVerbose>3)	cout << "Comparing with tube " << hit[ii][0];
	  Found = kFALSE;
	  if (neighborings.At(ineigh)==hit[ii][0]) {
	    if (fVerbose>3) cout << " ... yes!" << endl;
	    if (IsAdded[(int)hit[ii][0]][(int)hit[ii][1]] == kFALSE) {
	      ordMap[ncl].push_back(make_pair(hit[ii][0],hit[ii][1]));
	      jt = ordMap[ncl].begin();
	      IsAdded[(int)hit[ii][0]][(int)hit[ii][1]] = kTRUE;
	      if (fVerbose>2) cout << "==> Adding tube " << hit[ii][0] << "\t hit# " << hit[ii][1] << endl;
	      Found = kTRUE;
	    }
	    else {
	      if (ncl>0) {
		for (Int_t jj=0; jj<listmulti.GetSize(); jj++) {
		  if (hit[ii][0] == listmulti.At(jj)) {
		    Found = kFALSE;
		    break;
		  }
		  else 
		    Found = kTRUE;
		}
	      }
	      else Found = kTRUE;
	      if (fVerbose>3) cout << "Already added to the list" << endl;
	    }
	  }
	  else 
	    if (fVerbose>3) cout << " ... no! Going to the next one... " << endl;
	  if (Found==kTRUE)
	    break;
	}
      }
      IsChecked[itube][index] = kTRUE;
    }

    if (fVerbose>2) {
      cout << endl;
      cout << "==== Map after hit ordering [tubeID, hit#]" << endl;
      for (myiter it = ordMap[ncl].begin(); it!=ordMap[ncl].end(); it++)
	cout << "(" << it->first << ", " << it->second << ")  ";
      cout << endl;     cout << endl;
    }
    
    size[ncl] = ordMap[ncl].size();
    if (fVerbose>2) {
      cout << "Size of the ordered map: " << size[ncl] << endl;
      cout << "There are " << hits_left - size[ncl] << " hits left " << endl;
    }

    hits_left -= size[ncl];
    if (hits_left!=0)
      ncl++;
    
  } while (hits_left!=0);

  if (endloop-size[ncl]>0) {
    for (int i=0; i<endloop; i++) {
      if (IsAdded[(int)hit[i][0]][(int)hit[i][1]] == kFALSE && IsFake[(int)hit[i][0]][(int)hit[i][1]] == kFALSE) {
	fakelist.push_back(make_pair((int)hit[i][0], (int)hit[i][1]));
	IsFake[(int)hit[i][0]][(int)hit[i][1]] = kTRUE;
      }
    }
  }

  if (fVerbose>2) {
    if (ncl==0)
      cout << "There is 1 cluster" << endl;
    else 
      cout << "There are " << ncl+1 << " clusters" << endl;
    for (int icl=0; icl<ncl+1; icl++) {
      cout << "Cluster " << icl << " --> size: " << size[icl] << endl;
      for (myiter it = ordMap[icl].begin(); it!=ordMap[icl].end(); it++)
	cout << "(" << it->first << ", " << it->second << ")  ";
      cout << endl;     cout << endl;
    }
    cout << endl;
  }

  return ncl+1;

}



// ========== LinFitPerp =================================================
// function for the linear fit --- without errors 
// minimizing the perpendicular
Bool_t PndTrkPrototype::LinFitPerp(TMatrixT<Double_t> hit, Int_t size, Double_t p[2], Double_t &meand) {

  fIsHor = kFALSE;

  if (fVerbose>1)
    cout  << "=== PndTrkPrototype::LinFitPerp" << endl;

  // linear fit without errors
  // y(a,b) = a+b*z
  // a = ymed - b*zmed
  // b1 = -B+sqrt(B*B+1), b2 = -B-sqrt(B*B+1)
  // B = 0.5*(szz-syy)/szy
  // szy = sum(zi*yi)-n*zmed*ymed, szz = sum(zi^2-nzmed^2), syy = sum(yi^2)-n*ymed^2
  // zmed = sum(zi)/n, ymed = sum(yi)/n

  if ((size<5 && !fSmallButOk) || size>64) { // to be changed into 64
    for (int i=0; i<size; i++)
      IsAdded[(int)hit[i][0]][(int)hit[i][1]] = kFALSE;
    fNfailed++;
    if (fVerbose>2) 
      cout << "PndTrkPrototype::LinFitPerp --- less than 5 hits in the cluster " 
 	   << " --> no tracking!" << endl;
    return kFALSE;
  }

  // check if track is horizontal
  Int_t nhor = 0;
  for (Int_t i=0; i<size-1; i++)
    if (hit[i][3] == hit[i+1][3]) nhor++;
  if (nhor==size-1) {
    fIsHor = kTRUE;
    p[0] = hit[0][3];
    p[1] = 0;
    if (fVerbose>1) 
      cout << "Track parameters after prefit: a = " << p[0] << ", b = " << p[1] << endl;
    return kTRUE;
  }
  else {
    fIsHor = kFALSE;
    
    Double_t zmed=0., ymed=0., zy=0., zz=0., yy=0.;
    Double_t szz=0., syy=0., szy=0., B=0., a1=0., b1=0., a2=0., b2=0.; 
    Double_t sumdist1=0., sumdist2=0., meandist1=0., meandist2=0.;
    Double_t p1[2], p2[2];
    p1[0] = 0; p1[1] = 0; p2[0] = 0; p2[1] = 0; 
    
    for (Int_t i=0; i<size; i++) {
      zmed += hit[i][2];
      ymed += hit[i][3];
      zy   += hit[i][2]*hit[i][3];
      yy   += hit[i][3]*hit[i][3];
      zz   += hit[i][2]*hit[i][2];
    }
    zmed = zmed/size;
    ymed = ymed/size;
    szz  = zz - size*zmed*zmed;
    syy  = yy - size*ymed*ymed;
    szy  = zy - size*ymed*zmed;
    if (szy!=0) {
      B = 0.5*(szz-syy)/szy;
      b1 = -B + TMath::Sqrt(B*B+1);
      b2 = -B - TMath::Sqrt(B*B+1);
      a1 = ymed - b1*zmed;
      a2 = ymed - b2*zmed;
      p1[0] = a1; p1[1] = b1;
      p2[0] = a2; p2[1] = b2;
      for (Int_t i=0; i<size; i++) {
	sumdist1 += Point2Line(hit[i][2], hit[i][3], p1);
	sumdist2 += Point2Line(hit[i][2], hit[i][3], p2);
      }
      meandist1 = sumdist1/size;
      meandist2 = sumdist2/size;
      
      if (fVerbose>1) {
	cout << "=== Prefit (LinFitPerp): choice of parameters: " << endl;
	cout << "Option 1: a = " << p1[0] << ", b = " << p1[1] << ", dist " << meandist1 << endl;
	cout << "Option 2: a = " << p2[0] << ", b = " << p2[1] << ", dist " << meandist2 << endl;
      }
      
      if (meandist1 < meandist2) {
	p[0] = p1[0]; p[1] = p1[1];
	meand = meandist1; // mm
      }
      else {
	p[0] = p2[0]; p[1] = p2[1];
	meand = meandist2;
      }
      if (fVerbose>1) 
	cout << "Track parameters after prefit: a = " << p[0] << ", b = " << p[1] << endl;
      return kTRUE;
    }
    else {
      if (fVerbose>1)
	cout << "szy = 0!!!" << endl;
      return kFALSE;
    }
  }
}


// ========== MinuitFit ==================================================
// implementation of the minuit fit
Bool_t PndTrkPrototype::MinuitFit(Int_t flag, Int_t size, TMatrixT<Double_t> hit, Double_t *pinit, Double_t pfit[2], Double_t errpfit[2], Double_t &fval) {

  if (fVerbose>1)
    cout << "=== PndTrkPrototype::MinuitFit" << endl;

  //  cout << "Size " << size << " flag " << fSmallButOk << endl;

  if ((size==0) || (size<5 && !fSmallButOk) || (size>64)) { // to be changed into 64
    if (fVerbose>2) 
      cout << "PndTrkPrototype::LinFitPerp --- less than 5 hits in the cluster " 
 	   << " --> no tracking!" << endl;
    return kFALSE;
  }

  pfit[0] = -999; pfit[1] = -999;

  TMinuit minimizer(2);
  minimizer.SetPrintLevel(-1);

  if (fVerbose>2) {
    cout << "****** MINUIT *******" << endl;
    cout << "a SEED " << pinit[1] << endl;
    cout << "b SEED " << pinit[0] << endl;
    cout << "********************" << endl;
  }

  TMatrixT<Double_t> fitvector;
  size = SetupFitVector(size, hit, fitvector);

  if (flag==1)
    minimizer.SetFCN(fcnStraightLine);
  else if (flag==2 || flag==4)
    minimizer.SetFCN(fcnLinePoints_new);
  else if (flag==3)
    minimizer.SetFCN(fcnLinePoints_old);
  minimizer.DefineParameter(0, "p0", pinit[0], 0.1, 0., 0.); //0., 0.);
  minimizer.DefineParameter(1, "p1", pinit[1], 0.1, 0., 0.);
  minimizer.SetObjectFit(&fitvector);

  minimizer.SetMaxIterations(500);

  Int_t checkminuit = minimizer.Migrad();
  hstats->Fill(checkminuit);

  if (checkminuit==0) {
    Double_t chisquare;
    minimizer.GetParameter(0, pfit[0], errpfit[0]);
    minimizer.GetParameter(1, pfit[1], errpfit[1]);
    
    if (fVerbose>1) {
      cout << "RESULTS MINUIT FIT ****************** " << endl;
      cout << "a = " << pfit[1] << "\t err a = " << errpfit[1] << endl;
      cout << "b = " << pfit[0] << "\t err b = " << errpfit[0] << endl;
    }
    Double_t grad[2];
    minimizer.Eval(2,grad,fval,pfit,3);

    Int_t ndf = size - minimizer.GetNumFreePars();
    chisquare = fval*ndf;
    Double_t cumul = 1-TMath::Prob(chisquare, ndf);
    
    // filling histos
    if (flag==1)
      chi2r_prefit->Fill(fval);
    else if (flag==2) {
      chi2r_fit->Fill(fval);
      cumul_fit->Fill(cumul);
    }
    else if (flag==3) {
      chi2r_oldfit->Fill(fval);
      cumul_oldfit->Fill(cumul);
    }
    else if (flag==4) {
      chi2r_refit->Fill(fval);
      cumul_refit->Fill(cumul);
    }
    
    return kTRUE;
  }

  else
    return kFALSE;

}


// ========== SetupFitVector =============================================
Int_t PndTrkPrototype::SetupFitVector(Int_t size, TMatrixT<Double_t> hit, TMatrixT<Double_t> &fitvect) 
{
  fitvect.ResizeTo(size,8);
  int counter = 0;
  if (fVerbose>2) 
    cout << "=========== SETUP FITVECTOR =============" << endl;
  for (int i=0; i<size; i++) 
    {
      for (int j=0; j<8; j++)
	fitvect[counter][j] = hit[i][j];
      if (fVerbose>3) {
	cout << "Checking FitVector: pt = (" << fitvect[counter][0] << ", " 
 	     << fitvect[counter][1] << "), isochrone = (" << fitvect[counter][2]
 	     << " +/- " << fitvect[counter][3] << ") cm " << endl; 
      }
      counter++;
    }
  if (size != counter)
    fitvect.ResizeTo(counter,8);
  if (fVerbose>2) cout << "Counter: " << counter << endl;

  return counter;

}

// ========== FCN =======================================================
void fcnStraightLine(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  f = 0.;
  Double_t chis = 0.;

  TMatrixT<Double_t> *mtx = (TMatrixT<Double_t> *) gMinuit->GetObjectFit();
  Int_t hitcounter = mtx->GetNrows();

  for (int i=0; i<hitcounter; i++) {
 
    Double_t dist = TMath::Abs(par[1]*mtx[0][i][2] - mtx[0][i][3] + par[0]) / TMath::Sqrt(par[1]*par[1] + 1);
    if (mtx[0][i][5]!=0)
      chis += pow(dist-mtx[0][i][4], 2) / pow(mtx[0][i][5],2);
    else
      chis += pow(dist-mtx[0][i][4], 2) / 12.;
  }
  if ((hitcounter-2) != 0)
    f = chis/(hitcounter-2);
}



// ========== FCN =======================================================
void fcnLinePoints_new(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  f = 0.;
  Double_t chis = 0.;

  TMatrixT<Double_t> *mtx = (TMatrixT<Double_t> *) gMinuit->GetObjectFit();
  Int_t hitcounter = mtx->GetNrows();
  Double_t sigmatot2 = 0;

  for (int i=0; i<hitcounter; i++) {
    sigmatot2 = mtx[0][i][5]*mtx[0][i][5]*(1-par[1]*par[1])*(1-par[1]*par[1])/((1+par[1]*par[1])*(1+par[1]*par[1]));
    Double_t dist = TMath::Abs(mtx[0][i][3]-par[0]-par[1]*mtx[0][i][2])/TMath::Sqrt(1+par[1]*par[1]);
    if (sigmatot2) 
      chis += dist*dist / sigmatot2;
    else
      chis += dist*dist / 12.;
  }
  if ((hitcounter-2) != 0)
    f = chis/(hitcounter-2);
}

// ========== FCN =======================================================
void fcnLinePoints_old(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  f = 0.;
  Double_t chis = 0.;

  TMatrixT<Double_t> *mtx = (TMatrixT<Double_t> *) gMinuit->GetObjectFit();
  Int_t hitcounter = mtx->GetNrows();
  Double_t sigmatot2 = 0;

  for (int i=0; i<hitcounter; i++) {
    sigmatot2 = 2*par[1]*par[1]*mtx[0][i][5]*mtx[0][i][5]/((1+par[1]*par[1])*(1+par[1]*par[1]));
    Double_t dist = TMath::Abs(mtx[0][i][3]-par[0]-par[1]*mtx[0][i][2])/TMath::Sqrt(1+par[1]*par[1]);
    if (sigmatot2) 
      chis += dist*dist / sigmatot2;
    else
      chis += dist*dist / 12.;
  }
  if ((hitcounter-2) != 0)
    f = chis/(hitcounter-2);
}


// ========== IntersectionFinder =========================================
// function to find the coordinate of the hit on the isochrone
// the coordinates of this point are taken from the intersection between the 
// isochrone (circumference with radius obtained from drift time) 
// and the line orthogonal to the fitted track,
// passing through the center of the wire [zc, yc]
Bool_t PndTrkPrototype::IntersectionFinder(Int_t size, TMatrixT<Double_t> hit, Double_t* p, TMatrixT<Double_t> & intersection) 
{
  if (fVerbose>1)
    cout << "=== PndTrkPrototype::IntersectionFinder" << endl;

  for (int i=0; i<size; i++) {
    // 1. line passing through the wire center [zc, yc] 
    // orthogonal to the fitted track (y = p[0] + p[1]*z)
    // --> y = m*z + q
    // m = -1/p[1], q = -m*zc +zc
    
    Double_t m, q;
    if (p[1]==0) 
      return kFALSE;
    else {
      m = -1./p[1];
      q = -m*hit[i][2]+hit[i][3];
      
      // 2. intersection of the line in 1. with the isochrone
      // + and - refers to the 2 possible intersections
      // + --> (z1, y1)
      Double_t z1, y1, z2, y2;
      z1 = (-(m*(q - hit[i][3]) - hit[i][2]) + sqrt((m*(q - hit[i][3]) - hit[i][2])*(m*(q - hit[i][3]) - hit[i][2]) - (m*m + 1)*((q - hit[i][3])*(q - hit[i][3]) + hit[i][2]*hit[i][2] - hit[i][4]*hit[i][4]))) / (m*m + 1);
      y1 = m*z1 + q;
      // -  --> (z2,y2)
      z2 = (-(m*(q - hit[i][3]) - hit[i][2]) - sqrt((m*(q - hit[i][3]) - hit[i][2])*(m*(q - hit[i][3]) - hit[i][2]) - (m*m + 1)*((q - hit[i][3])*(q - hit[i][3]) + hit[i][2]*hit[i][2] - hit[i][4]*hit[i][4]))) / (m*m + 1);
      y2 = m*z2 + q;
      
      // 3. distance between point (z1, y1) and the fitted track (the same for (z2, y2))
      Double_t dist1 = TMath::Abs(p[1]*z1-y1+p[0])/TMath::Sqrt(p[1]*p[1]+1);
      Double_t dist2 = TMath::Abs(p[1]*z2-y2+p[0])/TMath::Sqrt(p[1]*p[1]+1);
      
      // 4. choice of the point --> hit on the isochrone (zh, yh);
      if (dist1 <= dist2) {
	intersection[i][2] = z1; 
	intersection[i][3] = y1;
      }
      else {
	intersection[i][2] = z2; 
	intersection[i][3] = y2;
      }
      for (int j=0; j<8; j++) 
	if (j!= 2 && j!=3)
	  intersection[i][j] = hit[i][j];
    }
  }
  return kTRUE;
}


// ========== Point2Line =================================================
Double_t PndTrkPrototype::Point2Line(Double_t x, Double_t y, Double_t *p) {
  return TMath::Abs(p[1]*x - y + p[0]) / TMath::Sqrt(p[1]*p[1] + 1);
}


// ========== CheckDistance ==============================================
// function to delete spurious hits from prefit on the basis of the hit
// distance from the prefit straight line
Bool_t PndTrkPrototype::CheckDistance(Int_t size, TMatrixT<Double_t> hit, Double_t *p, Int_t &fake, TMatrixT<Double_t> &newhit, vector<pair<int,int> > &fakelist) {

  if (fVerbose>1)
    cout << "=== PndTrkPrototype::CheckDistance " << endl;

  fake=0;
  Double_t distance;

  if (fVerbose>2) {
    cout << "SIZE " << size << endl;
    for (Int_t i=0; i<size; i++) 
      cout << hit[i][0] << "\t";
    cout << endl;
  }

  for (Int_t i=0; i<size; i++) {
    distance = Point2Line(hit[i][2], hit[i][3], p);
    hdist->Fill(distance);
    if (distance>0.8 && IsFake[(int)hit[i][0]][(int)hit[i][1]] == kFALSE) {
      fakelist.push_back(make_pair(hit[i][0], hit[i][1]));
      IsAdded[(int)hit[i][0]][(int)hit[i][1]] = kFALSE;
      IsChecked[(int)hit[i][0]][(int)hit[i][1]] = kFALSE;
      IsFake[(int)hit[i][0]][(int)hit[i][1]] = kTRUE;
      if (fVerbose>1)
	cout << "Tube " << hit[i][0] << " (hit " << hit[i][1] << ") -- distance from prefit: " << distance << " cm " << endl;
      for (Int_t nn=i; nn<size-1; nn++) {
	for (Int_t k=0; k<8; k++)
	  hit[nn][k] = hit[nn+1][k];
      }
      i--;
      size = size - 1;
      fake++;
      hit.ResizeTo(size,8);
    }
    if (fVerbose>3)
      cout << "After " << i+1 << " checks, the size is: " << size << endl;
  }

  if (fake==0)
    return kFALSE;

  if (fake>0) {
    newhit.ResizeTo(size,8);
    newhit = hit;
  }
  if (fVerbose > 1) {
    cout << "There were " << fake << " fake hits and the new size is " << size << endl;
    if (fake>0 && fVerbose>1) {
      cout << "The new hit matrix is: " << endl;
      for (int i=0; i<size; i++) 
	cout << newhit[i][0] << "\t";
      cout << endl;
    }
  }

  return kTRUE;
}


// ========== CheckResiduals =============================================
// function to delete spurious hits from prefit on the basis of the residuals
Bool_t PndTrkPrototype::CheckResiduals(Int_t flag, Int_t size, TMatrixT<Double_t> hit, Double_t *p, Int_t &fake, TMatrixT<Double_t> &newhit, vector<pair<int,int> > &fakelist) {

  if (fVerbose>1)
    cout << "=== PndTrkPrototype::CheckResiduals " << endl;

  fake=0;
  Double_t residual;
  Int_t init_size = size;
  TMatrixT<Double_t> init_hit;
  init_hit.ResizeTo(init_size,8);
  init_hit = hit;

  for (Int_t i=0; i<size; i++) {
    residual = Point2Line(hit[i][2], hit[i][3], p) - hit[i][4];
    if (flag==2)
      hresfit->Fill(residual*10); // mm 
    else if (flag==3) {
      hresrefit->Fill(residual*10); // mm
    }
    if (fVerbose>3)
      cout << "Hit number " << i << ": " << hit[i][2] << "\t res = " << residual << endl;
    if (TMath::Abs(residual)>0.2 && IsFake[(int)hit[i][0]][(int)hit[i][1]] == kFALSE) { // cm
      if (fVerbose>2)
	cout << "Residual for tube " << hit[i][0] << " is " << residual << " cm " << endl;
      fakelist.push_back(make_pair(hit[i][0], hit[i][1]));
      IsAdded[(int)hit[i][0]][(int)hit[i][1]] = kFALSE;
      IsChecked[(int)hit[i][0]][(int)hit[i][1]] = kFALSE;
      IsFake[(int)hit[i][0]][(int)hit[i][1]] = kTRUE;
      for (Int_t nn=i; nn<size-1; nn++) {
	for (Int_t k=0; k<8; k++)
	  hit[nn][k] = hit[nn+1][k];
      }
      i--;
      size = size - 1;
      fake++;
      hit.ResizeTo(size,8);
    }
  }
  if (fVerbose>2)
    cout << "After " << size << " checks, the size is: " << size << endl;

  if (fake==0)
    return kFALSE;

  if (fake>=2*init_size/3.) {
    if (fVerbose>2)
      cout << "%%%%%%%%%%%  # fake >= 2*(#hits)/3 [= " << 2*init_size/3. << "]" << endl;
    Double_t yl1 = init_hit[0][3];
    Double_t yl2 = init_hit[0][3];
    Int_t layer1 = 0, layer2 = 0;
    for (int i=0; i<init_size; i++) {
      if (init_hit[i][3] == yl1) layer1++;
      else {
	yl2 = init_hit[i][3];
	layer2++;
      }
    }
    if (fVerbose>2) 
      cout << "# hits in layer1: " << layer1 << ", # hits in layer2: " << layer2 << endl;
    if (layer1==16 && layer2==16) {
      Int_t count = 0;
      fakelist.clear();
      for (int i=0; i<init_size; i++) {
	if (init_hit[i][3]==yl1) {
	  newhit.ResizeTo(count+1,8);
	  for (int j=0; j<8; j++) 
	    newhit[count][j] = init_hit[i][j];
	  count++;
	}
	else if (init_hit[i][3]==yl2) {
	  fakelist.push_back(make_pair(init_hit[i][0], init_hit[i][1]));
	  IsAdded[(int)init_hit[i][0]][(int)init_hit[i][1]] = kFALSE;
	  IsChecked[(int)init_hit[i][0]][(int)init_hit[i][1]] = kFALSE;
	  IsFake[(int)init_hit[i][0]][(int)init_hit[i][1]] = kTRUE;
	}
      }
      size = layer1;
      fake = layer2;
    }
    else {
      if (size>0) {
	newhit.ResizeTo(size,8);
	newhit = hit;
      }
      else
	return kFALSE;
    }
  }
  else {
    if (size>0) {
      newhit.ResizeTo(size,8);
      newhit = hit;
    }
    else
      return kFALSE;
  }

  if (fVerbose > 2) {
    cout << "There were " << fake << " fake hits and the new size is " << size << endl;
    if (fVerbose>2 && fake>0 && size>0) {
      cout << "The new hit matrix is: " << endl;
      for (int i=0; i<size; i++) 
	cout << newhit[i][0] << "\t";
      cout << endl;
    }
  }

  if (fake>0 && init_size!=fake) 
    return kTRUE;
}


// ========== TrackIntersections =========================================
// function that checks which tubes the reconstructed tracks intersects
Int_t PndTrkPrototype::TrackIntersections(Int_t nhits, Double_t* p, TMatrixT<Double_t> hitOrig, TMatrixT<Double_t> hit, TMatrixT<Double_t>& newhypohit, Int_t &size, Int_t &k, Double_t &trklength, Double_t &esum) {

  if (fVerbose>1)
    cout << "=== PndTrkPrototype::TrackIntersections" << endl;

  Double_t R = 0.5;
  Double_t a = p[0];
  Double_t b = p[1];
  Double_t zc = 0., yc = 0.;
  Double_t A, B, C, Delta;
  Bool_t HasFired = kFALSE;
  Bool_t AlreadyAdded = kFALSE;
  Int_t wrong = 0;
  trklength = 0.;
  Double_t res = 0.;
  k = 0;
  PndSttTube* tube;

  // retrieving all tubes
  for(int itube = 1; itube < fSttTubeArray->GetEntriesFast(); itube++) {
    res = 0;
    AlreadyAdded = kFALSE;
    tube = (PndSttTube*) fSttTubeArray->At(itube);
    zc = tube->GetPosition().Z();
    yc = tube->GetPosition().Y();
    if (fVerbose>3) 
      cout << "Checking tube " << itube << endl;
    A = (1 + b*b);
    B = (zc - a*b + b*yc);
    C = zc*zc + a*a - 2*a*yc + yc*yc - R*R;
    Delta = B*B - A*C;
    if (fVerbose>3) cout << "DELTA " << Delta << endl;
    HasFired = kFALSE;
    if (Delta>=0) {
      if (fVerbose>2) cout << "Tube " << itube << " is hit by the recon track ";
      Double_t z1 = (B + TMath::Sqrt(Delta))/A;
      Double_t y1 = a + b*z1;
      Double_t z2 = (B - TMath::Sqrt(Delta))/A;
      Double_t y2 = a + b*z2;
      Double_t dist = TMath::Sqrt((z1-z2)*(z1-z2) + (y1-y2)*(y1-y2));
      trklength += dist;
      for (Int_t i=0; i<size; i++) {
	if (zc==hit[i][2] && yc==hit[i][3]) {
	  HasFired = kTRUE;
	  esum += hit[i][6];
	  if (fVerbose>2) cout << " AND it was hit by the real track" << endl;
	}
	if (HasFired==kTRUE) break;
      }
      if (HasFired == kFALSE) {
	if (fVerbose>2)
	  cout << " BUT it was NOT hit by the real track " << endl;
	wrong++;
      }
      HasFired = kFALSE;
      for (int jj=0; jj<fNSttHits; jj++) {
	if(zc==hitOrig[jj][2] && yc==hitOrig[jj][3]) {
	  for (int jjj=0; jjj<k; jjj++) {
	    if (newhypohit[jjj][0] == hitOrig[jj][0]) {
	      AlreadyAdded = kTRUE;
	      break;
	    }
	  }
	  if (AlreadyAdded==kFALSE && IsUsed[(int)hitOrig[jj][1]]==kFALSE) {
	    newhypohit.ResizeTo(k+1,8);
	    for (int kk=0; kk<8; kk++) 
	      newhypohit[k][kk] = hitOrig[jj][kk];
	    if (fVerbose>2)
	      cout << "Adding to newhypo hit # " << hitOrig[jj][1] << ", tube " << hitOrig[jj][0] << endl;
	    k++;
	  }
	}
      }
    }
  }

  for (int i=0; i<size; i++) {
    for (int j=0; j<k; j++) {
      if (hit[i][0] == newhypohit[j][0]) {
	IsAdded[(int)hit[i][0]][(int)hit[i][1]] = kTRUE;
	break;
      }
      else 
	IsAdded[(int)hit[i][0]][(int)hit[i][1]] = kFALSE;
    }
  }

  if (fVerbose>2)
    cout << "There are " << wrong << " additional hits!!!" << endl;

  return wrong;
}


// =======================================================================
// =========================  Drawing functions  =========================
// =======================================================================
// ========== DrawGeometry ===============================================
void PndTrkPrototype::DrawGeometry()
{
 
  h = new TH2F("zy","zy", 100, 470, 670, 100, -100, 100);
  c->cd();
  h->Draw();
  h->SetStats(kFALSE);

  TArc *arc = NULL;
  TMarker *mrk = NULL;
  PndSttTube *tube = NULL;

  // draw all the tubes
  for(int itube = 1; itube < fSttTubeArray->GetEntriesFast(); itube++) {
    tube = (PndSttTube*) fSttTubeArray->At(itube);
    arc = new TArc(tube->GetPosition().Z()*10, tube->GetPosition().Y()*10, 0.5*10);
    arc->SetFillStyle(0);
    arc->Draw("SAME");
  }
  c->Update();
  c->Modified();
}


// ========== DrawIntersectionPoints =====================================
void PndTrkPrototype::DrawMCPoint(TMatrixT<Double_t> points, Int_t size, Double_t m, Double_t q)
{
  TMarker *mrk;
  c->cd();
  
  for (Int_t i=0; i<size; i++) {
    mrk = new TMarker(points[i][0]*10, points[i][1]*10, 3);
    mrk->SetMarkerSize(1.5);
    mrk->Draw("SAME");
  }
  Double_t pts[2][2];
  pts[0][0] = 490;
  pts[1][0] = 670;

  TLine *l = new TLine(pts[0][0], m*pts[0][0]+q, pts[1][0], m*pts[1][0]+q);
  l->Draw("same");

  c->Update();
  c->Modified();
}


// ========== DrawCluster ================================================
void PndTrkPrototype::DrawCluster(Int_t iter, TMatrixT<Double_t> hit, Int_t size)
{
  TMarker *m;
  c->cd();
  for (Int_t i=0; i<size; i++) {
    m = new TMarker(hit[i][2]*10, hit[i][3]*10, 7);
    if (iter==0)
      m->SetMarkerColor(37);
    else if (iter==1)
      m->SetMarkerColor(kCyan);
    m->Draw("SAME");
  }
  c->Update();
  c->Modified();
}


// ========== DrawIsochrone ==============================================
void PndTrkPrototype::DrawIsochrone(TMatrixT<double> hit, Int_t size, Int_t iter, Int_t ncluster)
{
  TArc *isochrone;
  TMarker *m;
  c->cd();

  if (fVerbose>2) 
    cout << "==== DrawIsochrone --- iter " << iter << ", cluster # " << ncluster << endl;
  for (Int_t i=0; i<size; i++) {
    if (fVerbose>3)
      cout << "Drawing tube " << hit[i][0] << ", hit " << hit[i][1] << endl;
    isochrone = new TArc(hit[i][2]*10, hit[i][3]*10, hit[i][4]*10);
    m = new TMarker(hit[i][2]*10, hit[i][3]*10, 7);
    if (iter==0){
      isochrone->SetLineColor(37);
      m->SetMarkerColor(37);
    }
    else if (iter==-1) {
      isochrone->SetLineColor(kBlack);
      m->SetMarkerColor(kBlack);
    }
    else {
      if (ncluster==0) {
	isochrone->SetLineColor(kMagenta-7);
	m->SetMarkerColor(kMagenta-7);
      }
      else {
	isochrone->SetLineColor(kCyan+1);
	m->SetMarkerColor(kCyan+1);
      }
    }
    isochrone->SetFillStyle(0);
    isochrone->Draw("SAME");
    m->Draw("SAME");
  }
  c->Update();
  c->Modified();
}


// ========== DrawIntersectionPoints =====================================
void PndTrkPrototype::DrawIntersections(TMatrixT<Double_t> hit, Int_t size, Int_t iter, Int_t ncluster)
{
  TMarker *m;
  c->cd();
  for (Int_t i=0; i<size; i++) {
    m = new TMarker(hit[i][2]*10, hit[i][3]*10, 5);
    if (iter==0)
      m->SetMarkerColor(kBlue);
    else if (iter==1) {
      if (ncluster==0)
	m->SetMarkerColor(kMagenta);
      else 
	m->SetMarkerColor(kBlue);
    }
    m->Draw("SAME");
  }

  c->Update();
  c->Modified();
}


// ========== DrawFit ====================================================
void PndTrkPrototype::DrawFit(Int_t step, Int_t iter, Double_t *p, Int_t ncluster)
{
  c->cd();

  if (fVerbose>2)
    if (step==1)
      cout << "=== PndTrkPrototype::DrawPrefit --> a = " << p[0] << ", b = " << p[1] << endl;
    else if (step==2)
      cout << "=== PndTrkPrototype::DrawMinuitFit --> a = " << p[0] << ", b = " << p[1] << endl;
    else if (step==3)
      cout << "=== PndTrkPrototype::DrawMinuitFit2 --> a = " << p[0] << ", b = " << p[1] << endl;

  TLine *l;
  Double_t pts[2][2];
  pts[0][0] = 490;
  pts[1][0] = 670;
  pts[0][1] = (p[0] + p[1]*pts[0][0]/10)*10;
  pts[1][1] = (p[0] + p[1]*pts[1][0]/10)*10;
  l = new TLine(pts[0][0], pts[0][1], pts[1][0], pts[1][1]);
  switch (step) {
  case 1: 
    if (iter==0) 
      l->SetLineColor(kRed); 
    else if (iter==1)
      l->SetLineColor(kCyan);
    break;
  case 2:
    if (iter==0)
      l->SetLineColor(kBlue); 
    else if (iter==1)
      l->SetLineColor(kCyan+3); 
    break;
  case 3:
    if (iter==0)
      l->SetLineColor(kGreen); 
    else if (iter==1) {
      l->SetLineColor(kGreen+2);
    }
    else if (iter==2) {
      if (ncluster==0)
	l->SetLineColor(kMagenta);
      else 
	l->SetLineColor(kBlue);

    }
    break;
  }
  l->Draw("SAME");

  c->Update();
  c->Modified();
}


// ========== WriteHistograms ============================================
void PndTrkPrototype::WriteHistograms() {

  if (fVerbose>0) {
    cout << endl;
    cout << "=== PndTrkPrototype::WriteHistograms " << endl;
    cout << endl;
    cout << "SUMMARY VALUES" << endl;
    cout << "Reco " << nReco << endl;
    cout << "Reco Less " << nRecoLess << "\t reco more " << nRecoMore << endl;
    cout << nReco << " reconstructed tracks out of " << nTotalMC << " MC tracks with STT hits " << endl;
    cout << "--> Efficiency = " << (double)100*nReco/nTotalMC << "%" << endl;
    cout << endl;
  }

  TFile* file = FairRootManager::Instance()->GetOutFile();
  file->cd();
  file->mkdir("PndTrkPrototype");
  file->cd("PndTrkPrototype");
  hdist->Write();
  hresfit->GetXaxis()->SetTitle("Residuals (mm)");
  hresfit->Write();
  hresrefit->GetXaxis()->SetTitle("Residuals (mm)");
  hresrefit->Write();
  hlastresid->GetXaxis()->SetTitle("Residuals (mm)");
  hlastresid->Write();
  chi2r_prefit->GetXaxis()->SetTitle("#chi^{2}_{r}");
  chi2r_prefit->Write();
  chi2r_fit->GetXaxis()->SetTitle("#chi^{2}_{r}");
  chi2r_fit->Write();
  chi2r_refit->GetXaxis()->SetTitle("#chi^{2}_{r}");
  chi2r_refit->Write();
  chi2r_oldfit->GetXaxis()->SetTitle("#chi^{2}_{r}");
  chi2r_oldfit->Write();
  cumul_fit->GetXaxis()->SetTitle("Cumulative distribution");
  cumul_fit->Write();
  cumul_refit->GetXaxis()->SetTitle("Cumulative distribution");
  cumul_refit->Write();
  cumul_oldfit->GetXaxis()->SetTitle("Cumulative distribution");
  cumul_oldfit->Write();
  hstats->Write();
  hmeanres->GetXaxis()->SetTitle("Mean residual/track (mm)" );
  hmeanres->Write();
  htrklength->Write();
  hde->Write();
  hdedx->Write();
  htubemult->Write();
  hgoodd->Write();
}


Double_t PndTrkPrototype::GetTimeInfo(Double_t radius)
{
  TF1 *frt = new TF1("frt", "[0]+[1]*x+[2]*pow(x,2)+[3]*pow(x,3)+[4]*pow(x,4)", 0, 500);

  if (fFlag==1)
    frt->SetParameters(0.0034801, 0.0534479, -6.32247e-5, -1.11069e-6, 4.31749e-9);
  else if (fFlag==2)
    frt->SetParameters(0.0034801, 0.0534479, -6.32247e-5, -1.11069e-6, 4.31749e-9);

  Double_t time = frt->GetX(radius*10);
  delete frt;
  return time;
}



ClassImp(PndTrkPrototype)