//  #include "glpk.h"
#include "PndTrkTracking.h"
#include "PndTrkBoundaryParStraws.h"
#include "PndTrkChi2Fits.h"
#include "PndTrkComparisonMCtruth.h"
#include "PndTrkSttConformalFilling.h"
// #include "PndTrkGlpkFits.h"
#include "PndTrkLegendreFits.h"
#include "PndTrkCleanup.h"
#include "PndTrkCTFindTrackInXY.h"
#include "PndTrkCTGeometryCalculations.h"
#include "PndTrkMergeSort.h"
#include "PndTrkPlotMacros.h"
#include "PndTrkPrintouts.h"
#include "PndTrkSttAdjacencies.h"

#include "PndSttHit.h"
#include "PndSciTHit.h"
#include "PndSttPoint.h"
#include "PndSttTrack.h"
#include "PndSttPoint.h"
#include "PndSttHelixHit.h"
#include "PndSttSingleStraw.h"
#include "PndSttTube.h"
#include "PndSttMapCreator.h"

#include "PndSdsHit.h"
#include "PndSdsMCPoint.h"

#include "PndTrackCand.h"
#include "PndTrackCandHit.h"
#include "PndTrack.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 <iostream>
#include <cmath>




// definition of constants used in this class

// integer constants
#define MAXHITSINCELL		20
// the following MAXHITSINFIT cannot be too large because the fit
// crashes 'silently' for too much memory consumption in the character arrays
// or takes too long time;


// MAXHITSINFIT has to be 12 when using the GLPK fitting procedure; there is no limitation
// with the other type of fits;
//#define MAXHITSINFIT		12
#define MAXHITSINFIT		50
#define MAXMVDMCPOINTS		2000
#define MAXSKEWHITSINFIT	8
#define MINIMUMHITSPERTRACK	3
#define MINOUTERHITSPERTRACK	5
#define TIMEOUT			60



// floating point constants

#define APOTEMAMAXINNERPARSTRAW	23.246827
#define APOTEMAMAXSKEWSTRAW	31.517569 // delimitation of the skew area
#define APOTEMAMINOUTERPARSTRAW	31.863369
#define APOTEMAMINSKEWSTRAW	23.246827 // delimitation of the skew area
#define BFIELD			2.  // in Tesla
#define CVEL			2.99792  //  velocity of light
#define DELTAnR			2.   //range of nR in TrkAssociatedParallelHitsToHelixquater
#define DIAMETERSTRAWTUBE	1.
#define DIMENSIONSCITIL		2.85 // cm
#define ERRORPIXEL		0.02611
#define ERRORSTRIP		0.02611
#define ERRORSQPIXEL		0.00068175
#define ERRORSQSTRIP		0.00068175
#define PI			3.141592654
#define PMAX			100.
#define RSTRAWDETECTORMAX	40.73 // maximum radius of the Stt detector in  cm
#define RSTRAWDETECTORMIN	16.119 // minimum radius of the Stt detector in  cm
#define STRAWRADIUS		0.5
#define STRAWRESOLUTION		0.015
#define STRAW_SKEW_INCLINATION	3.
#define STTDRIFTVEL		0.0025	//   in cm/nsec
#define VERTICALGAP		4. // (cm) gap between Left and Right sections of detector.



using namespace std;


// -----   Default constructor   -------------------------------------------
PndTrkTracking::PndTrkTracking() : FairTask("Tracking") { 
 fPersistence = kTRUE;
 istampa = 0;
 iplotta = false;
 doMcComparison = false;
 fYesClean = false;
 fYesCleanMvd = true;
 fYesSciTil = false ;
 fMvdAloneTracking = false;
 fNevents_to_plot = 10;


 Initialization_ClassVariables();
 sprintf(fSttBranch,"STTHit");
 sprintf(fMvdPixelBranch,"MVDHitsPixel");
 sprintf(fMvdStripBranch,"MVDHitsStrip");
}
// -------------------------------------------------------------------------

PndTrkTracking::PndTrkTracking(Int_t verbose) : FairTask("Tracking") { 
 fPersistence = kTRUE;
 istampa = verbose;
 iplotta = false;
 doMcComparison = false;
 fYesClean = false;
 fYesCleanMvd = true;
 fYesSciTil = false ;
 fMvdAloneTracking = false;

 Initialization_ClassVariables();
 sprintf(fSttBranch,"STTHit");
 sprintf(fMvdPixelBranch,"MVDHitsPixel");
 sprintf(fMvdStripBranch,"MVDHitsStrip");
}
// -------------------------------------------------------------------------

PndTrkTracking::PndTrkTracking(int istamp, bool  iplot, bool imc)
		: FairTask("Tracking") { 
 fPersistence = kTRUE;
 istampa = istamp;
 iplotta = iplot;
 doMcComparison = imc;
 fYesClean = false;
 fYesCleanMvd = true;
 fYesSciTil = false ;
 fMvdAloneTracking = false;
 Initialization_ClassVariables();
 sprintf(fSttBranch,"STTHit");
 sprintf(fMvdPixelBranch,"MVDHitsPixel");

 sprintf(fMvdStripBranch,"MVDHitsStrip");
}


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

PndTrkTracking::PndTrkTracking(int istamp, bool  iplot, bool imc, bool doSciTil)
		: FairTask("Tracking") { 
 fPersistence = kTRUE;
 istampa = istamp;
 iplotta = iplot;
 doMcComparison = imc;
 fYesClean = false;
 fYesCleanMvd = true;
 fYesSciTil = doSciTil ;
 fMvdAloneTracking = false;
 Initialization_ClassVariables();
 sprintf(fSttBranch,"STTHit");
 sprintf(fMvdPixelBranch,"MVDHitsPixel");

  sprintf(fMvdStripBranch,"MVDHitsStrip");
}
// -----   Destructor   ----------------------------------------------------
PndTrkTracking::~PndTrkTracking() {}
// -------------------------------------------------------------------------

//--------------- begin  PndTrkTracking::Initialization_ClassVariables

void PndTrkTracking::Initialization_ClassVariables()
{
// this is only for initializing the Class Variables.

	size_t len;


// booleans :
	len = sizeof(SingleHitListStt);
	memset (SingleHitListStt,true,len);

	len = sizeof(fInclusionListSciTil);
	memset (fInclusionListSciTil,true,len);

	len = sizeof(fInclusionListStt);
	memset (fInclusionListStt,true,len);

	len = sizeof(finMvdTrackCandPixel);
	memset (finMvdTrackCandPixel,false,len);

	len = sizeof(finMvdTrackCandStrip);
	memset (finMvdTrackCandStrip,false,len);

	len = sizeof(fTypeConf);
	memset (fTypeConf,false,len);

// char :


	len = sizeof(fSttBranch);
	memset (fSttBranch,0,len);

	len = sizeof(fMvdPixelBranch);
	memset (fMvdPixelBranch,0,len);

	len = sizeof(fMvdStripBranch);
	memset (fMvdStripBranch,0,len);

//  Short_t :

	fnMCTracks=0;
	fnSciTilHits=0;

	len = sizeof(fListMvdPixelHitsinTrack);
	memset (fListMvdPixelHitsinTrack,0,len);

	len = sizeof(fListMvdStripHitsinTrack);
	memset (fListMvdStripHitsinTrack,0,len);

	len = sizeof(fListSciTilHitsinTrack);
	memset (fListSciTilHitsinTrack,0,len);

	len = sizeof(fListSttParHits);
	memset (fListSttParHits,0,len);

	len = sizeof(fListSttParHitsinTrack);
	memset (fListSttParHitsinTrack,0,len);

	len = sizeof(fListSttSkewHitsinTrack);
	memset (fListSttSkewHitsinTrack,0,len);

	len = sizeof(fListSttSkewHits);
	memset (fListSttSkewHits,0,len);

	len = sizeof(fListSttSkewHitsinTrackSolution);
	memset (fListSttSkewHitsinTrackSolution,0,len);

	len = sizeof(fListTrackCandHit);
	memset (fListTrackCandHit,0,len);

	len = sizeof(fnMvdPixelHitsinTrack);
	memset (fnMvdPixelHitsinTrack,0,len);

	len = sizeof(fnMvdStripHitsinTrack);
	memset (fnMvdStripHitsinTrack,0,len);

	len = sizeof(fnTrackCandHit);
	memset (fnTrackCandHit,0,len);

	len = sizeof(fnSciTilHitsinTrack);
	memset (fnSciTilHitsinTrack,0,len);

	len = sizeof(fnSttParHitsinTrack);
	memset (fnSttParHitsinTrack,0,len);

	len = sizeof(fnSttSkewHitsinTrack);
	memset (fnSttSkewHitsinTrack,0,len);

	fnMvdDSPixelHitNotTrackCand=0;
	fnMvdDSStripHitNotTrackCand=0;
	fnMvdPixelHit=0;
	fnMvdStripHit=0;
	fnMvdTrackCand=0;
	fnMvdUSPixelHitNotTrackCand=0;
	fnMvdUSStripHitNotTrackCand=0;

	len = sizeof(fnHitMvdTrackCand);
	memset (fnHitMvdTrackCand,0,len);

	len = sizeof(fListHitMvdTrackCand);
	memset (fListHitMvdTrackCand,0,len);

	len = sizeof(fListHitTypeMvdTrackCand);
	memset (fListHitTypeMvdTrackCand,0,len);

	len = sizeof(fListMvdDSPixelHitNotTrackCand);
	memset (fListMvdDSPixelHitNotTrackCand,0,len);

	len = sizeof(fListMvdUSPixelHitNotTrackCand);
	memset (fListMvdUSPixelHitNotTrackCand,0,len);

	len = sizeof(fListMvdDSStripHitNotTrackCand);
	memset (fListMvdDSStripHitNotTrackCand,0,len);

	len = sizeof(fListMvdUSStripHitNotTrackCand);
	memset (fListMvdUSStripHitNotTrackCand,0,len);

	len = sizeof(fListTrackCandHitType);
	memset (fListTrackCandHitType,0,len);

	len = sizeof(fnParContiguous);
	memset (fnParContiguous,0,len);

	len = sizeof(fListParContiguous);
	memset (fListParContiguous,0,len);

//  int :
	fNevents_to_plot = 10;


//  Double_t :

	fFimin=0.;
	SEMILENGTH_STRAIGHT=75.;
	ZCENTER_STRAIGHT=35.;

	len = sizeof(fALFA);
	memset (fALFA,0,len);

	len = sizeof(fBETA);
	memset (fBETA,0,len);

	len = sizeof(fGAMMA);
	memset (fGAMMA,0,len);

	len = sizeof(fCxMC);
	memset (fCxMC,0,len);

	len = sizeof(fCyMC);
	memset (fCyMC,0,len);

	len = sizeof(fR_MC);
	memset (fR_MC,0,len);

	len = sizeof(fMCtruthTrkInfo);
	memset (fMCtruthTrkInfo,0,len);

	len = sizeof(fMCSkewAloneX);
	memset (fMCSkewAloneX,0,len);

	len = sizeof(fMCSkewAloneY);
	memset (fMCSkewAloneY,0,len);

	len = sizeof(fradiaConf);
	memset (fradiaConf,0,len);

	len = sizeof(fOx);
	memset (fOx,0,len);

	len = sizeof(fOy);
	memset (fOy,0,len);

	len = sizeof(fR);
	memset (fR,0,len);

	len = sizeof(frefindexMvdPixel);
	memset (frefindexMvdPixel,0,len);

	len = sizeof(fsigmaXMvdPixel);
	memset (fsigmaXMvdPixel,0,len);

	len = sizeof(fsigmaYMvdPixel);
	memset (fsigmaYMvdPixel,0,len);

	len = sizeof(fsigmaZMvdPixel);
	memset (fsigmaZMvdPixel,0,len);

	len = sizeof(fXMvdPixel);
	memset (fXMvdPixel,0,len);

	len = sizeof(fYMvdPixel);
	memset (fYMvdPixel,0,len);

	len = sizeof(fZMvdPixel);
	memset (fZMvdPixel,0,len);

	len = sizeof(fXMvdStrip);
	memset (fXMvdStrip,0,len);

	len = sizeof(fYMvdStrip);
	memset (fYMvdStrip,0,len);

	len = sizeof(fZMvdStrip);
	memset (fZMvdStrip,0,len);

	len = sizeof(fsigmaXMvdStrip);
	memset (fsigmaXMvdStrip,0,len);

	len = sizeof(fsigmaYMvdStrip);
	memset (fsigmaYMvdStrip,0,len);

	len = sizeof(fsigmaZMvdStrip);
	memset (fsigmaZMvdStrip,0,len);

	len = sizeof(frefindexMvdStrip);
	memset (frefindexMvdStrip,0,len);

	len = sizeof(fposizSciTil);
	memset (fposizSciTil,0,len);



//  pointers :

	HANDLE=NULL;
	HANDLE2=NULL;

	hdeltaRPixel=NULL;
	hdeltaRStrip=NULL;
	hdeltaRPixel2=NULL;
	hdeltaRStrip2=NULL;
	fMCTrackArray=NULL;
	fSttTubeArray=NULL;
	fSttPointArray=NULL;
	fSttHitArray=NULL;
	fSttTrackArray=NULL;
	fSttTrackCandArray=NULL;
	fMvdPixelHitArray=NULL;
	fMvdStripHitArray=NULL;
	fMvdTrackCandArray=NULL;
	fSciTHitArray=NULL;
	fMvdMCPointArray=NULL;
	fSttMvdPndTrackCandArray=NULL;
	fSttMvdPndTrackArray=NULL;
	fSttParameters=NULL;

}

//--------------- end  PndTrkTracking::Initialization_ClassVariables


// -----   Public method Init   --------------------------------------------

InitStatus PndTrkTracking::Init() {


 IVOLTE=-1;

 SEMILENGTH_STRAIGHT = 75.;
 ZCENTER_STRAIGHT = 35.;


 if(iplotta){
  hdeltaRPixel = new TH1F("hdeltaRPixel", "distance MC Pixel point from trajectory in XY plane", 100, -1, 1);
  hdeltaRStrip = new TH1F("hdeltaRStrip", "distance MC Strip point from trajectory in XY plane", 100, -1, 1);
  hdeltaRPixel2 = new TH1F(
  "hdeltaRPixel2", "distance MC point from trajectory in XY plane (Pixels)", 100, -10, 10);
  hdeltaRStrip2 = new TH1F(
  "hdeltaRStrip2", "distance MC point from trajectory in XY plane (Strips)", 100, -10,10);
}


//  --------------------------- opening files for special purposes


if(doMcComparison >=1 ){
//---- apertura file con info su Found tracce su cui si fa Helix fit dopo
   HANDLE2 = fopen("info_da_PndTrackFinderReal.txt","w");

//  ---- open filehandle per statistica sugli hits etc.
   HANDLE = fopen("statistiche.txt","w");
//  ---------------


}  //  end of if(istampa >=1)


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



 // Get RootManager
 FairRootManager* ioman = FairRootManager::Instance();
// ioman = FairRootManager::Instance();
 if ( ! ioman ) {
    cout << "-E- PndTrkTracking::Init: "
	 << "RootManager not instantiated, return!" << endl;
    return kFATAL;
 }
//  -----   maps of STT tubes
 // CHECK added 
 PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters);
 fSttTubeArray = mapper->FillTubeArray();
 //----------------------------------------------------  end map

// load the array indicating if a straw is exetrnal of not;
// true -->  it is external; false --> it is internal;
// remember that the numbering of the STT Straws starts at 1;
	PndTrkBoundaryParStraws BoundaryParStraws ;
	fExternal_Straws[0] = false;  // just to be super-safe;
	for(int i=1;i<= NUMBER_STRAWS;i++){
		fExternal_Straws[i] = BoundaryParStraws.Set(i);
	}
 //----------------------------------------------------


//    get   the MCTrack  array

 fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
 if ( ! fMCTrackArray)
 {
	cout << "-E- PndTrkTracking::Init: No MCTrack array, return!"
	   << endl;
	return kERROR;
 }

//  -------------------------   get the SciTil hits
 if(fYesSciTil) {
	fSciTHitArray = (TClonesArray*) ioman->GetObject("SciTHit");
 } else {
	fSciTHitArray = NULL;
 }
//---------------------------

//  -------------------------   get the SciTil MC Points
 if(fYesSciTil && doMcComparison) {
	fSciTPointArray = (TClonesArray*) ioman->GetObject("SciTPoint");
 } else {
	fSciTPointArray = NULL;
 }
//---------------------------
 // Get input array   questi sono i MC point di STT
 fSttPointArray = (TClonesArray*) ioman->GetObject("STTPoint");
 if ( ! fSttPointArray ) {
	cout << "-W- PndSttHelixHitProducer::Init: "
	 << "No STTPoint array, return!" << endl;
	return kERROR;
 }

 // Get input array   hit di STT after digi
 fSttHitArray = (TClonesArray*) ioman->GetObject(fSttBranch);
//  fSttHitArray = (TClonesArray*) ioman->GetObject("STTHit");
 if ( ! fSttHitArray ) {
	cout << "-W- PndTrkTracking::Init: "
	 << "No STTHit array, return!" << endl;
	return kERROR;
 }

//  -------------------------   get the Mvd hits
 fMvdPixelHitArray = (TClonesArray*) ioman->GetObject(fMvdPixelBranch);
//  fMvdPixelHitArray = (TClonesArray*) ioman->GetObject("MVDHitsPixel");
 if ( !fMvdPixelHitArray){
	cout << "-W- PndTrkTracking::Init: " << "No MVD Pixel hitArray!" <<endl;

//	return kERROR;
 }
 fMvdStripHitArray = (TClonesArray*) ioman->GetObject(fMvdStripBranch);
//  fMvdStripHitArray = (TClonesArray*) ioman->GetObject("MVDHitsStrip");

 if ( !fMvdStripHitArray){
	cout << "-W- PndTrkTracking::Init: " << "No MVD Strip hitArray!" <<endl;
//	return kERROR;
 }

//  -------------------------   get the Mvd track candidates

 if ( fMvdAloneTracking ){
	fMvdTrackCandArray = (TClonesArray*) ioman->GetObject("MVDRiemannTrackCand");
	if ( !fMvdTrackCandArray){
		cout << "-W- PndTrkTracking::Init: " << "No MVD TrackCand Array, return!" <<endl;
		return kERROR;
	}
 }
 cout << "-I- PndTrkTracking: Initialization successfull" << endl;

//  -------------------------   get the Mvd MC points

 fMvdMCPointArray = (TClonesArray*) ioman->GetObject("MVDPoint");
 if ( !fMvdMCPointArray){
	cout << "-W- PndTrkTracking::Init: " << "No MVD MC Point Array!" <<endl;
//	return kERROR;
 }
 cout << "-I- PndTrkTracking: Initialization successfull" << endl;


//--------------------------------  output TClonesArrays ------------------------

 // Create and register output array for PndTrackCand of Stt+Mvd combined

 fSttMvdPndTrackCandArray = new TClonesArray("PndTrackCand");
 ioman->Register("SttMvdTrackCand","SttMvd",fSttMvdPndTrackCandArray, fPersistence);


 // Create and register output array for PndTrack of Stt+Mvd combined

 fSttMvdPndTrackArray = new TClonesArray("PndTrack");
 ioman->Register("SttMvdTrack","SttMvd",fSttMvdPndTrackArray, fPersistence);

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


//   calculate the boundaries of the Box in Conformal Space, see Gianluigi logbook on pag. 210-211
 Short_t i;
 Double_t
	A,
	r1,
	r2;

 fradiaConf[0] = 1./RSTRAWDETECTORMAX;
 r1 = RSTRAWDETECTORMIN;
 A = (RSTRAWDETECTORMAX - r1)/NRDIVCONFORMAL;
 if ( NRDIVCONFORMAL > 1 ) {
	for(i = 1; i< NRDIVCONFORMAL ; i++){
		r2 = r1 + A;
		fradiaConf[NRDIVCONFORMAL-i] = 1./r2;
		r1=r2;
	}
 }

 // load the adjacencies table for the Stt tubes;
 PndTrkSttAdjacencies  Adjacent;
 Adjacent.CalculateAdjacentStt(
	NUMBER_STRAWS,
	fSttTubeArray,
	fnParContiguous,  // output; number of contiguous straws (axial Stt);
	fListParContiguous  // output list (axial Stt);
	);


 return kSUCCESS;

}

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




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




void PndTrkTracking::WriteHistograms(){

  TFile* file = FairRootManager::Instance()->GetOutFile();
  file->cd();
  file->mkdir("PndTrkTracking");
  file->cd("PndTrkTracking");
  hdeltaRPixel->Write();
  hdeltaRStrip->Write();
  hdeltaRPixel2->Write();
  hdeltaRStrip2->Write();
  delete hdeltaRPixel;
  delete hdeltaRStrip;
  delete hdeltaRPixel2;
  delete hdeltaRStrip2;

}

// -----   Public method Exec   --------------------------------------------
// -----   Public method Exec   --------------------------------------------
void PndTrkTracking::Exec(Option_t* opt) {


 bool
	flag,
	intersect,
	outcome;

//----------------------------------------------
 bool
	tGoodSkewFit[MAXTRACKSPEREVENT],
	tkeepit[MAXTRACKSPEREVENT],
	tMvdhits[MAXTRACKSPEREVENT],
	tstatus[MAXTRACKSPEREVENT],
	tSttSZfit[MAXTRACKSPEREVENT];

	Vec <bool>
	GoodSkewFit(tGoodSkewFit,MAXTRACKSPEREVENT,"GoodSkewFit"),
	keepit(tkeepit ,MAXTRACKSPEREVENT,"keepit"),
	Mvdhits(tMvdhits, MAXTRACKSPEREVENT,"Mvdhits"),
	status(tstatus, MAXTRACKSPEREVENT,"status"),
	SttSZfit(tSttSZfit, MAXTRACKSPEREVENT,"SttSZfit");

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

memset(tMvdhits,false,sizeof(tMvdhits));





 Short_t
	tnHitsInMCTrack[MAXTRACKSPEREVENT],
	tnMCParalAlone[MAXTRACKSPEREVENT],
	tnMCSkewAlone[MAXTRACKSPEREVENT],
	tnParalCommon[MAXTRACKSPEREVENT],
	tnSkewCommon[MAXTRACKSPEREVENT],
	tnSkewHitsInMCTrack[MAXTRACKSPEREVENT],
	tnSpuriParinTrack[MAXTRACKSPEREVENT],
	tnSpuriSkewinTrack[MAXTRACKSPEREVENT],
	//  given a Hit number it gives its radial box number
	tRConformalIndex[MAXSTTHITS],
	//  given a Hit number it gives its azimuthal box number
	tFiConformalIndex[MAXSTTHITS],
	ttempore[MAXSTTHITS],
	TemporarySkewList[2*MAXSTTHITS][2],
	BigList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK],
	// nBoxConformal,  first index -> radial divisions,
	// 2nd index -> azimuthal divisions; n. of hits falling in this cell.
	tnBoxConformal[NRDIVCONFORMAL*NFIDIVCONFORMAL],
	tParalCommonList[MAXTRACKSPEREVENT*MAXSTTHITSINTRACK],
	tParSpuriList[MAXTRACKSPEREVENT*MAXSTTHITSINTRACK],
	tSkewCommonList[MAXTRACKSPEREVENT*MAXSTTHITSINTRACK],
	tSkewSpuriList[MAXTRACKSPEREVENT*MAXSTTHITSINTRACK],
	tHitsinBoxConformal[MAXHITSINCELL*NRDIVCONFORMAL*NFIDIVCONFORMAL],
	tCharge[MAXTRACKSPEREVENT],
	tdaTrackFoundaTrackMC[MAXTRACKSPEREVENT],
	tresultFitSZagain[MAXTRACKSPEREVENT],
	tstatusflag[MAXTRACKSPEREVENT];

 Vec <Short_t>
	nHitsInMCTrack(tnHitsInMCTrack, MAXTRACKSPEREVENT,"nHitsInMCTrack"),
	nMCParalAlone(tnMCParalAlone, MAXTRACKSPEREVENT,"nMCParalAlone"),
	nMCSkewAlone(tnMCSkewAlone, MAXTRACKSPEREVENT,"nMCSkewAlone"),
	nParalCommon(tnParalCommon,MAXTRACKSPEREVENT,"nParalCommon"),
	nSkewCommon(tnSkewCommon,MAXTRACKSPEREVENT,"nSkewCommon"),
	nSkewHitsInMCTrack(tnSkewHitsInMCTrack, MAXTRACKSPEREVENT,"nSkewHitsInMCTrack"),
	nSpuriParinTrack(tnSpuriParinTrack, MAXTRACKSPEREVENT,"nSpuriParinTrack"),
	nSpuriSkewinTrack(tnSpuriSkewinTrack, MAXTRACKSPEREVENT,"nSpuriSkewinTrack"),
	//  given a Hit number it gives its radial box number
	RConformalIndex(tRConformalIndex, MAXSTTHITS,"RConformalIndex"),
	//  given a Hit number it gives its azimuthal box number
	FiConformalIndex(tFiConformalIndex, MAXSTTHITS,"FiConformalIndex"),
	tempore(ttempore, MAXSTTHITS,"tempore"),
	// nBoxConformal,  first index -> radial divisions,
	// 2nd index -> azimuthal divisions; n. of hits falling in this cell.
	nBoxConformal(tnBoxConformal, NRDIVCONFORMAL*NFIDIVCONFORMAL,"nBoxConformal"),
	ParalCommonList(tParalCommonList, MAXTRACKSPEREVENT*MAXSTTHITSINTRACK,"ParalCommonList"),
	ParSpuriList(tParSpuriList, MAXTRACKSPEREVENT*MAXSTTHITSINTRACK,"ParSpuriList"),
	SkewCommonList(tSkewCommonList, MAXTRACKSPEREVENT*MAXSTTHITSINTRACK,"SkewCommonList"),
	SkewSpuriList(tSkewSpuriList, MAXTRACKSPEREVENT*MAXSTTHITSINTRACK,"SkewSpuriList"),
	HitsinBoxConformal(tHitsinBoxConformal, MAXHITSINCELL*NRDIVCONFORMAL*NFIDIVCONFORMAL,"HitsinBoxConformal"),
	Charge(tCharge,MAXTRACKSPEREVENT,"Charge"),
	daTrackFoundaTrackMC(tdaTrackFoundaTrackMC,MAXTRACKSPEREVENT,"daTrackFoundaTrackMC"),
	resultFitSZagain(tresultFitSZagain,MAXTRACKSPEREVENT,"resultFitSZagain"),
	statusflag(tstatusflag,MAXTRACKSPEREVENT,"statusflag");

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




 Short_t
	nalone,
	ncand,
	nhitsinfit,
	Nint,
	NNN,
	nRemainingCandidates,
	npixelhitsintrack,
	nstriphitsintrack,
	nTotalCandidates,
	nXYZhits,
	i,
	iParHit,
	ipunto,
	j,
	k,
	kall,
	l
	;


 Int_t
	iaccept,
	len,
	nSttHit,
	nSttMCPoint,
	nSttParHit,
	nSttSkewHit,
	nSttTrackCand
	;

 Int_t	nrounds0,
	nrounds1,
	tListHits[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK];

 Vec <Int_t> ListHits(tListHits,MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK,"ListHits");

 Double_t
	tS[2*MAXSTTHITS],
	tZ[2*MAXSTTHITS];
 Vec <Double_t>
	S(tS,2*MAXSTTHITS,"S"),
	Z(tZ,2*MAXSTTHITS,"Z");

 Double_t
	Distance,
	Fi,
	Ptras,
	ddd,
	delta,
	dis,
	dista,
	dista0,
	dista1,
	emme,
	gap,
	highqualitycut,
	Pxini,
	px,
	Pyini,
	py,
	Pzini,
	qop,
	x,
	y,
//
	AloneX[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK],
	AloneY[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK],
	DriftRadiusbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+
		MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK],// all skew hits have double
	ErrorchosenPixel[MAXMVDPIXELHITS],
	ErrorchosenStrip[MAXMVDSTRIPHITS],
	ErrorchosenSkew[MAXSTTHITS],
	ErrorDriftRadiusbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+
				MAXMVDSTRIPHITSINTRACK+2],// solutions
	FI0[MAXTRACKSPEREVENT],
	Fi_final_helix_referenceframe[MAXTRACKSPEREVENT],
	Fi_initial_helix_referenceframe[MAXTRACKSPEREVENT],
	Fi_low_limit[MAXTRACKSPEREVENT],
	Fi_up_limit[MAXTRACKSPEREVENT],
	KAPPA[MAXTRACKSPEREVENT],
	Posiz1[3],
	primoangolo[MAXTRACKSPEREVENT],
	s[2],
	Sbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK], // multiplication by 2 in the
	SchosenPixel[MAXTRACKSPEREVENT][MAXMVDPIXELHITS],
	SchosenStrip[MAXTRACKSPEREVENT][MAXMVDSTRIPHITS],
	SchosenSkew[MAXTRACKSPEREVENT][MAXSTTHITS], // NO multiplication by 2 here because for the
					// skew hits only one
					// solution is selected.
	Start[3],

	TemporaryS[2*MAXSTTHITS],
	TemporaryZ[2*MAXSTTHITS],
	TemporaryZDrift[2*MAXSTTHITS],
	TemporaryZErrorafterTilt[2*MAXSTTHITS],

	temporeZErrorafterTilt[MAXSTTHITSINTRACK],
	temporeS[MAXSTTHITSINTRACK],
	temporeZ[MAXSTTHITSINTRACK],
	temporeZDrift[MAXSTTHITSINTRACK],
	tmpErrorZDrift[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK],
	tmpS[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK],
	tmpZ[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK],
	tmpZDrift[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK],
	trajectory_vertex[2],
	Trajectory_Start[MAXTRACKSPEREVENT][2],
	ultimoangolo[MAXTRACKSPEREVENT],
	versor[2],
	WDX[MAXSTTHITS],
	WDY[MAXSTTHITS],
	WDZ[MAXSTTHITS],
	X[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK],
	XintersectionList[2],
	Y[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK],
	YintersectionList[2],
	z[2],
	ZEDbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK], // rather improbable chance that
	zdrift[2],
	zerror[2],
	ZDrift[2*MAXSTTHITS],
	ZErrorafterTilt[2*MAXSTTHITS],
	zeta0,
	zeta1,
//
	info[MAXSTTHITS][7],
	infoparalConformal[MAXSTTHITS][5],
	Sfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK],
	XY[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK][2],
	ZchosenPixel[MAXTRACKSPEREVENT][MAXMVDPIXELHITS],
	ZchosenStrip[MAXTRACKSPEREVENT][MAXMVDSTRIPHITS],
	ZchosenSkew[MAXTRACKSPEREVENT][MAXSTTHITS],
	Zfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK],
	ZDriftfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK],
	ZErrorafterTiltfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK];


//--------------------
	Double_t
	tDriftRadius[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK],
	tErrorDriftRadius[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK],
	tZED[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK];

	Vec <Double_t>
	DriftRadius(tDriftRadius,MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK,"DriftRadius"),
	ErrorDriftRadius(tErrorDriftRadius,MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK,"ErrorDriftRadius"),
	ZED(tZED,MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK
		+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK,"ZED");
//--------------------------------------

 TVector3
	ErrMomentum,
	ErrPosition,
	Momentum,
	Position;

 FairMCPoint *puntator;

 PndMCTrack* pMCtr;

 PndSttHit *pSttHit;

 PndSttTube *pSttTube;

 PndSdsHit
	*pMvdPixelHit,
	*pMvdStripHit;

 PndTrkCleanup Cleaner;

 PndTrkCTGeometryCalculations GeomCalculator;

 PndTrkPrintouts fPrint;


 // the class with all the fits. This is used for the SZ fit.
 bool YesGLPKfitSZ = false;
// PndTrkGlpkFits fit;  YesGLPKfitSZ=true;
// PndTrkLegendreFits fit;
 PndTrkChi2Fits fit;

//  reset the TClones Arrays of the PndTrackCand and PndTrack; it is necessary
//  to do this for every event at the very beginning of the Exec (those TClones Arrays
//  have been registered  in the Init and they are automatically writted each event).

 fSttMvdPndTrackCandArray->Delete();
 fSttMvdPndTrackArray->Delete();

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

 IVOLTE++;

 if(istampa>0)
	cout<<endl<<"Entering in PndTrkTrack : evt (starting from 0)  n. "<<IVOLTE<<endl;

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

// -------------------------------------  fetch info from MVD

 if( fMvdPixelHitArray ) {
	fnMvdPixelHit = fMvdPixelHitArray->GetEntriesFast();
 } else {
 	fnMvdPixelHit = 0;
 }

 if( fMvdStripHitArray ) {
	fnMvdStripHit = fMvdStripHitArray->GetEntriesFast();
 } else {
	fnMvdStripHit = 0;
 }

 if(fnMvdPixelHit>MAXMVDPIXELHITS){
	cout<<"from PndTrkTracking, fnMvdPixelHit is > maximum allowed ("
	<<MAXMVDPIXELHITS<<") and therefore is set to "
	<<MAXMVDPIXELHITS<<endl;
	fnMvdPixelHit=MAXMVDPIXELHITS;
 }
 if(fnMvdStripHit>MAXMVDSTRIPHITS){
	cout<<"from PndTrkTracking, fnMvdStripHit is > maximum allowed ("
	<<MAXMVDSTRIPHITS<<") and therefore is set to "
	<<MAXMVDSTRIPHITS<<endl;
	fnMvdStripHit=MAXMVDSTRIPHITS;
 }

// ---------------------------------------------  extract info from HITS Pixel MVD

 for( i= 0; i< fnMvdPixelHit; i++){
	pMvdPixelHit = (PndSdsHit *) fMvdPixelHitArray->At(i);
	TVector3 temp = pMvdPixelHit->GetPosition();
	fXMvdPixel[i] = temp.X();
	fYMvdPixel[i] = temp.Y();
	fZMvdPixel[i] = temp.Z();
	fsigmaXMvdPixel[i] = pMvdPixelHit->GetDx();
	fsigmaYMvdPixel[i] = pMvdPixelHit->GetDy();
	fsigmaZMvdPixel[i] = pMvdPixelHit->GetDz();
	frefindexMvdPixel[i] = pMvdPixelHit->GetRefIndex();
 }

// -------------------------------------------  extract info from HITS Strip MVD

 for( i= 0; i< fnMvdStripHit; i++){
	pMvdStripHit = (PndSdsHit *) fMvdStripHitArray->At(i);
	TVector3 temp = pMvdStripHit->GetPosition();
	fXMvdStrip[i] = temp.X();
	fYMvdStrip[i] = temp.Y();
	fZMvdStrip[i] = temp.Z();
	fsigmaXMvdStrip[i] = pMvdStripHit->GetDx();
	fsigmaYMvdStrip[i] = pMvdStripHit->GetDy();
	fsigmaZMvdStrip[i] = pMvdStripHit->GetDz();
	frefindexMvdStrip[i] = pMvdStripHit->GetRefIndex();
 }


// ---------------   printout of Mvd Hits;
 if(istampa>=1)
	fPrint.stampaMvdHits(
		fMvdPixelBranch,
		fMvdStripBranch,
		fnMvdPixelHit,
		fnMvdStripHit,
		frefindexMvdPixel,
		frefindexMvdStrip,
		fsigmaXMvdPixel,
		fsigmaXMvdStrip,
		fsigmaYMvdPixel,
		fsigmaYMvdStrip,
		fsigmaZMvdPixel,
		fsigmaZMvdStrip,
		fXMvdPixel,
		fXMvdStrip,
		fYMvdPixel,
		fYMvdStrip,
		fZMvdPixel,
		fZMvdStrip
			);


// ------------------------------------------ get info from trackcand  of MVD

 if( fMvdAloneTracking ){
	fnMvdTrackCand = fMvdTrackCandArray->GetEntriesFast();
	if (fnMvdTrackCand> MAXMVDTRACKSPEREVENT) {
		cout<<"from PndTrkTracking  :  N. of MvdTrackCand = "<<
		fnMvdTrackCand<<" and it is > MAXMVDTRACKSPEREVENT (="<<MAXMVDTRACKSPEREVENT
		<<"),  therefore it is set to "<<MAXMVDTRACKSPEREVENT<<endl;
		fnMvdTrackCand= MAXMVDTRACKSPEREVENT;
		ExtractInfoFromMvdTrackCand();
	}

 } else {
	fnMvdTrackCand=0;
 }





//---------------------------------------------   fetching the STT  MC points
 nSttMCPoint = fSttPointArray->GetEntriesFast();
 if (nSttMCPoint ==0){
	if(istampa>0) cout<<"warning from PndTrkTracking  :  N. of Stt MC points = 0"<<endl;
//	return;
 } else  if( nSttMCPoint>MAXSTTHITS){
	cout<<"warning from PndTrkTracking  :  N. of Stt MC points = "<<nSttMCPoint
	<<" and it is > MAXSTTHITS ("<<MAXSTTHITS
	     <<"), therefore consider only the first "<<MAXSTTHITS<<" hits"<<endl<<endl;
	nSttMCPoint= MAXSTTHITS;
 }


//--------------------------------------------------------   fetching the STT  hits

 nSttHit = fSttHitArray->GetEntriesFast();
 if (nSttHit ==0){
	if (istampa >= 1) cout<<"warning from PndTrkTracking, evt "<<IVOLTE<<"  :  N. di Stt Hits = 0, return!"<<endl<<endl;
	return;
 } else if (nSttHit> MAXSTTHITS) {
	cout<<"warning from PndTrkTracking, evt "<<IVOLTE<<"  :  N. of Stt Hits = "<<nSttHit
	     <<" and it is > MAXSTTHITS (="<<MAXSTTHITS
	     <<"), therefore consider only the first "<<MAXSTTHITS<<" hits"<<endl<<endl;
	nSttHit= MAXSTTHITS;
 }

 if (istampa >= 1) {
	cout<<"from PndTrkTracking  :  n. total # Hits in STT  : "<<nSttHit<<endl;
 }

 nSttParHit=0; 
 nSttSkewHit=0;

 for( i= 0; i< nSttHit; i++){
	pSttHit = (PndSttHit *) fSttHitArray->At(i);
	// right way to extract the corrisponding MC point.
	ipunto= pSttHit->GetRefIndex();
//	tubeID = pSttHit->GetTubeID();
	// fTubeID[i] = STT tubeID correspondint to the hit number i ;
	fTubeID[i] = pSttHit->GetTubeID();
	pSttTube = (PndSttTube *) fSttTubeArray->At(fTubeID[i]);
	// drift radius
	Double_t dradius = pSttHit->GetIsochrone();
	// wire direction
	TVector3 wiredirection = pSttTube->GetWireDirection();

	if(wiredirection.Z() >=0.) {
		WDX[i] = wiredirection.X();
		WDY[i] = wiredirection.Y();
		WDZ[i] = wiredirection.Z();
	} else {
		WDX[i] = -wiredirection.X();
		WDY[i] = -wiredirection.Y();
		WDZ[i] = -wiredirection.Z();
	}
	info[i][0]= pSttTube->GetPosition().X();
	info[i][1]= pSttTube->GetPosition().Y();
	info[i][2]= pSttTube->GetPosition().Z();
	info[i][3]= dradius;
	info[i][4]= pSttTube->GetHalfLength();

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


	if(ipunto>=0) {
		puntator = (FairMCPoint*) fSttPointArray->At(ipunto);
		info[i][6]= puntator->GetTrackID();
	} else {
		info[i][6]= -10.;
	}

	if( fabs( WDX[i] )< 0.00001 && fabs( WDY[i] )< 0.00001 ){
		info[i][5]= 1.;
		fListSttParHits[nSttParHit]=i;
		nSttParHit++;
	} else {
		info[i][5]= 99.;// to signal that it is a skew straw.
		fListSttSkewHits[nSttSkewHit]=i;
		nSttSkewHit++;
	}

 //  printout of the Stt hits;
 if (istampa >= 1 ) fPrint.stampaSttHits2(i,ipunto,dradius,WDX,WDY,WDZ,puntator,pSttTube,fTubeID[i]);

  }  //   end of for( i= 0; i< nSttHit; i++)

//   reordering the list of parallel hits ( fListSttParHits) by decreasing spatial radius;
//   first the outermost then the innermost. This is necessary because later the search
//   must starts from the outer hits. 

 Initial_SttParHits_DecreasingR_Ordering( info, fListSttParHits,nSttParHit );


//	fill the inclusion list for Stt, include only first hit for those straws with
//	multiple hits.

 MakeInclusionListStt(nSttHit, info);

//-----------------------------------  end of exclusion of straws with multiple hits



//-------------------------------------------- fetch the SciTil hits
 fnSciTilHits = 0;
 if( fSciTHitArray != NULL){
	// number SciTil hits/event
	fnSciTilHits = fSciTHitArray->GetEntriesFast();
	if(fnSciTilHits>MAXSCITILHITS){
		cout<<"from PndTrkTracking  :  N. of SciTil Hits = "<<fnSciTilHits
	     <<" and it is > MAXSCITILHITS (="<<MAXSCITILHITS
	     <<"), therefore consider only the first "<<MAXSCITILHITS <<" hits"<<endl<<endl;
		fnSciTilHits= MAXSCITILHITS;
	}
 }  // end of if( fSciTHitArray != NULL)

 // it is important that fSciTilMaxNumber and OriginalSciTilList have a scope
 // extending in all   exec  method; that's why they are stated here, out
 // of the if( fSciTHitArray != NULL).
 Short_t fSciTilMaxNumber;
 if( fnSciTilHits > 0) fSciTilMaxNumber = fnSciTilHits ; else fSciTilMaxNumber = 1;

 Short_t
	nHitsInSciTile[fSciTilMaxNumber],
	OriginalSciTilList[fSciTilMaxNumber][fSciTilMaxNumber];

//--------   initialization (to 0) of nHitsInSciTile array;
 memset (nHitsInSciTile,0,sizeof(nHitsInSciTile));
//---

 if(istampa>0) cout<<"da PndTrkTracking, event "<<IVOLTE<<", "<<fnSciTilHits
		<<" SciTil hits presenti inizialmente."<<endl;
	for(j=0; j<fnSciTilHits; j++){
	}

   if( fnSciTilHits >0 ){
	// OriginalSciTilList is the list of original SciTil hits (not purged yet)
	// present in a given SciTil tile :
	// OriginalSciTilList[nacceptedhit][*];


	PndSciTHit *pPndSciTHit;
	TVector3  posiz;
	// the first SciTil hit; this cannot be duplicate hit by definition.
	// The Sci Tiles are numbered here according to the numbering
	// of the (first) SciTil Hit inside the Sci Tile.
	pPndSciTHit = (PndSciTHit*) fSciTHitArray->At(0);
	posiz = pPndSciTHit->GetPosition();
	if(istampa>0)cout<<"da PndTrkTracking SciTil n. "<<0<<
	" non purgato, Xpos "<<
			posiz.X()<<", Ypos "<<
			posiz.Y()<<", Zpos "<<posiz.Z()<<endl;

	fposizSciTil[0][0]=posiz.X();
	fposizSciTil[0][1]=posiz.Y();
	fposizSciTil[0][2]=posiz.Z();
	fpSciTilx[0]=posiz.X();
	fpSciTily[0]=posiz.Y();
	fpSciTilz[0]=posiz.Z();
	OriginalSciTilList[0][0]=0;
	nHitsInSciTile[0]=1;
	iaccept=1;
	// the other SciTil hits; purge them if they are duplicate.
	for(j=1; j<fnSciTilHits; j++){
		pPndSciTHit = (PndSciTHit*) fSciTHitArray->At(j);
		posiz = pPndSciTHit->GetPosition();
		if(istampa>0)cout<<"da PndTrkTracking SciTil n. "<<j<<" non purgato, Xpos "
			<<posiz.X()<<", Ypos "<<posiz.Y()<<", Zpos "<<posiz.Z()<<endl;

	// purging the duplicate SciTil hits.
	    flag = true;
	    for(k=0; k<iaccept; k++){
		if(
			(fabs(posiz.X() - fposizSciTil[k][0])< 1.e-20) &&
			(fabs(posiz.Y() - fposizSciTil[k][1])< 1.e-20) &&
			(fabs(posiz.Z() - fposizSciTil[k][2])< 1.e-20)
		  ){
			flag=false;
			OriginalSciTilList[k][nHitsInSciTile[k]]= j;
			nHitsInSciTile[k]++;
			break;
		}  // end of if((fabs(posiz.X() - old...
	    } // end of for(k=0; k<iaccept; k++)
	    if(flag){
		fposizSciTil[iaccept][0]=posiz.X();
		fposizSciTil[iaccept][1]=posiz.Y();
		fposizSciTil[iaccept][2]=posiz.Z();
		fpSciTilx[iaccept]=posiz.X();
		fpSciTily[iaccept]=posiz.Y();
		fpSciTilz[iaccept]=posiz.Z();
		OriginalSciTilList[iaccept][0]= j;
		nHitsInSciTile[iaccept]=1;
		iaccept++;
	    }
	 }  // end of for(j=0; j<fnSciTilHits; j++)
	 fnSciTilHits=iaccept;

	 //  set the inclusion list of the SciTils to true.

	for(i=0;i<fnSciTilHits;i++){
		fInclusionListSciTil[i]=true;
	}


//-----------stampe.
if(istampa>0){
  cout<<"da PndTrkTracking, dopo purga di SciTil; n. hits = "<<fnSciTilHits<<endl;
  for(j=0; j<fnSciTilHits; j++){
	cout<<"da PndTrkTracking SciTil Xpos "<<fposizSciTil[j][0]<<", Ypos "<<
	fposizSciTil[j][1]<<", Zpos "<<fposizSciTil[j][2]<<endl;
  }
}
//---------- fine stampe.



   }  // end of if( fnSciTilHits >0 )


//-----------------------------------end fetching SciTil hits.

//-------------------
//-------------------
//-------------------
//-------------------
//-------------------
//-------------------
//-------------------  start the combined Mvd-Stt  PR
//-------------------
//-------------------
//-------------------
//-------------------
//-------------------


 //  initialization of the (assumed) starting point of the tracks, at the origin;
 //  needs to be modified later;

 trajectory_vertex[0]=trajectory_vertex[1]=0.;
 len = sizeof(Trajectory_Start);
 memset (Trajectory_Start,0,len);
//--------------


 PndTrkSttConformalFilling fill;

 fill.FromXYtoConformal(
	trajectory_vertex,
	info,
	fListSttParHits,
	nSttParHit,
	infoparalConformal,
	STRAWRADIUS
		);

 fill.BoxConformalFilling(
	FiConformalIndex,
	HitsinBoxConformal,
	fInclusionListStt,
	infoparalConformal,
	fListSttParHits,
	MAXHITSINCELL,
	nBoxConformal,
	NFIDIVCONFORMAL,
	nSttParHit,
	NRDIVCONFORMAL,
	fradiaConf,
	RConformalIndex
			);


 nSttTrackCand=0;    // # tracks found

//-----
 // class that finds the track (Stt hits only) in XY projection
 PndTrkCTFindTrackInXY SttTrackXYFinder;

 //  struct  necessary to pass all the parametrs to the PndTrkCTFindTrackInXY::FindTrackInXYProjection
 //  method. Since these parameters are > 60, cint does NOT accept to pass them in the usual
 //  way (parameters in the calling sequence) to PndTrkCTFindTrackInXY::FindTrackInXYProjection.
 FindTrackInXYProjection_InputData input;

 // loading those elements of the struct common to all the candidate tracks.
 input.apotemamaxskewstraw = APOTEMAMAXSKEWSTRAW;
 input.deltanr = DELTAnR;
 input.dimensionscitil = DIMENSIONSCITIL;
 input.FiConformalIndex = tFiConformalIndex;
 input.HitsinBoxConf = tHitsinBoxConformal;
 input.InclusionListStt = fInclusionListStt;
 input.InclusionListSciTil = fInclusionListSciTil;
 input.info = info;
 input.infoparalConformal = infoparalConformal;
 input.ListSttParHits = fListSttParHits;
 input.maxhitsinfit = MAXHITSINFIT;
 input.maxscitilhitsintrack = MAXSCITILHITSINTRACK;
 input.maxstthits = MAXSTTHITS;
 input.maxstthitsintrack = MAXSTTHITSINTRACK;
 input.minimumhitspertrack = MINIMUMHITSPERTRACK;
 input.minouterhitspertrack = MINOUTERHITSPERTRACK;
 input.nBoxConf = tnBoxConformal;
 input.nfidivconformal = NFIDIVCONFORMAL;
 input.nrdivconformal = NRDIVCONFORMAL;
 input.nSciTilHits = fnSciTilHits;
 input.nsttparhit = nSttParHit;
 input.posizSciT = fposizSciTil;
 input.radiaConf = fradiaConf;
 input.RConformalIndex = tRConformalIndex;
 input.rstrawdetectormax = RSTRAWDETECTORMAX;
 input.rstrawdetectormin = RSTRAWDETECTORMIN;
 input.strawradius = STRAWRADIUS;
 input.trajectory_vertex = trajectory_vertex;
 input.YesSciTil = fYesSciTil;

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





//----- loop over the SciTil hits first;


 //  U and V only for the parallel Stt hits.

 Double_t
	U[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK],
	V[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK];

 //  variable usedd in the plotting in the Legiandre plot;
 input.icounter=0;

 if(fYesSciTil) {
    for(i=0; i<fnSciTilHits ; i++) {
	if( nSttTrackCand >= MAXTRACKSPEREVENT ){
		cout<<"from PndTrkTracking :  # n. Tracks found so far = "
		<<nSttTrackCand<<" and it is >= MAXTRACKSPEREVENT ( = "
		 <<MAXTRACKSPEREVENT<<"); exiting from SciTil loop.\n";
		break ;
	}

	// inputs for the FindTrackInXYProjection class;
	input.iHit = -i-1;	// seed hit; it is negative for SciTil Hits.
	input.nRcell = -1;  // fR cell of the seed hit negative because SciTil hit is outside
			    // of the Stt system.;

	Fi =  atan2(fposizSciTil[i][1],fposizSciTil[i][0]) ;
//-----------------
	if ( Fi < 0. ) Fi += 2.*PI;
	input.nFicell =  (Short_t) (0.5*NFIDIVCONFORMAL*Fi/PI);
	if(input.nFicell > NFIDIVCONFORMAL ) {
		input.nFicell = NFIDIVCONFORMAL;
	} else if (input.nFicell<0) {
		input.nFicell = 0;
	}
	input.posizSciTilx = fposizSciTil[i][0];
	input.posizSciTily = fposizSciTil[i][1];

	// outputs from the FindTrackInXYProjection class are stored here;
	input.ALFA = &fALFA[nSttTrackCand];
	input.BETA = &fBETA[nSttTrackCand];
	input.Charge = &Charge[nSttTrackCand];
	input.Fi_final_helix_referenceframe = &Fi_final_helix_referenceframe[nSttTrackCand];
	input.Fi_initial_helix_referenceframe = &Fi_initial_helix_referenceframe[nSttTrackCand];
	input.Fi_low_limit = &Fi_low_limit[nSttTrackCand];
	input.Fi_up_limit = &Fi_up_limit[nSttTrackCand];
	input.GAMMA = &fGAMMA[nSttTrackCand];
	input.ListHitsinTrack = &fListSttParHitsinTrack[nSttTrackCand][0];
	input.ListSciTilHitsinTrack = &fListSciTilHitsinTrack[nSttTrackCand][0];
	input.nHitsinTrack = &fnSttParHitsinTrack[nSttTrackCand];
	input.nSciTilHitsinTrack = &fnSciTilHitsinTrack[nSttTrackCand];
	input.Oxx = &fOx[nSttTrackCand];
	input.Oyy = &fOy[nSttTrackCand];
	input.Rr = &fR[nSttTrackCand];
	input.S_SciTilHitsinTrack = &fS_SciTilHitsinTrack[nSttTrackCand][0];
	input.TypeConf = &fTypeConf[nSttTrackCand];
	input.U = &U[nSttTrackCand][0];
	input.V = &V[nSttTrackCand][0];

	input.icounter++;  // this is the plot number;
	outcome = SttTrackXYFinder.FindTrackInXYProjection(&input);

	if(!outcome){
		continue;
	}
	for(j=0; j<fnSttParHitsinTrack[nSttTrackCand]; j++){
	  fInclusionListStt[fListSttParHitsinTrack[nSttTrackCand][j]] = false;
	}
	keepit[nSttTrackCand]=true;
	nSttTrackCand++;

    } // end of  for(i=0; i<fnSciTilHits ; i++)


 }  // end of  if(fYesSciTil)


//----- end use the SciTil hits first.



//----- loop over the parallel hits

//   begins the first iteration with more severe cuts on the # hits in track candidate

int iconta=0;
// for(iParHit=0; iParHit<nSttParHit + 1 -  MINIMUMHITSPERTRACK ; iParHit++) {
 for(iParHit=0; iParHit<nSttParHit ; iParHit++) {
	if( ! fInclusionListStt[fListSttParHits[iParHit]] )  continue;
	if( !fExternal_Straws[fTubeID[fListSttParHits[iParHit]]] ) continue; // only seeds at the external boundary of the STT

	if( nSttTrackCand >= MAXTRACKSPEREVENT) {
		cout<<"from PndTrkTracking :  # n. Tracks found so far = "
		<<nSttTrackCand<<" and it is >= MAXTRACKSPEREVENT ( = "
		 <<MAXTRACKSPEREVENT<<"); exiting from || hit loop.\n";
		break;
	}
iconta++;
	// inputs for the FindTrackInXYProjection class;
	input.iHit = iParHit;// seed hit in the PARALLEL number scheme; it is negative for SciTil Hits.
	input.nFicell = FiConformalIndex[fListSttParHits[iParHit]];  // Fi cell of the seed hit;
	input.nRcell = RConformalIndex[fListSttParHits[iParHit]];  // fR cell of the seed hit;

	// outputs from the FindTrackInXYProjection class are stored here;
	input.ALFA = &fALFA[nSttTrackCand];
	input.BETA = &fBETA[nSttTrackCand];
	input.Charge = &Charge[nSttTrackCand];
	input.Fi_final_helix_referenceframe = &Fi_final_helix_referenceframe[nSttTrackCand];
	input.Fi_initial_helix_referenceframe = &Fi_initial_helix_referenceframe[nSttTrackCand];
	input.Fi_low_limit = &Fi_low_limit[nSttTrackCand];
	input.Fi_up_limit = &Fi_up_limit[nSttTrackCand];
	input.GAMMA = &fGAMMA[nSttTrackCand];
	input.ListHitsinTrack = &fListSttParHitsinTrack[nSttTrackCand][0];
	input.ListSciTilHitsinTrack = &fListSciTilHitsinTrack[nSttTrackCand][0];
	input.nHitsinTrack = &fnSttParHitsinTrack[nSttTrackCand];
	input.nSciTilHitsinTrack = &fnSciTilHitsinTrack[nSttTrackCand];
	input.Oxx = &fOx[nSttTrackCand];
	input.Oyy = &fOy[nSttTrackCand];
	input.Rr = &fR[nSttTrackCand];
	input.S_SciTilHitsinTrack = &fS_SciTilHitsinTrack[nSttTrackCand][0];
	input.TypeConf = &fTypeConf[nSttTrackCand];
	input.U = &U[nSttTrackCand][0];
	input.V = &V[nSttTrackCand][0];

	input.icounter++;  // this is the plot number;
	outcome = SttTrackXYFinder.FindTrackInXYProjection(&input);

	if(!outcome)  continue;


// --------  here the track and its hits were found, filling the Inclusion list

 for(j=0; j<fnSttParHitsinTrack[nSttTrackCand]; j++){
   fInclusionListStt[fListSttParHitsinTrack[nSttTrackCand][j]] = false;
 }

 keepit[nSttTrackCand]=true;
 nSttTrackCand++;

  }      // end  of   for(iParHit=0; iParHit<nSttParHit+1-MINIMUMHITSPERTRACK; iParHit++)



 for(i=0; i<nSttTrackCand;i++){

	//  set FI0 now; FI0 may be changed later in case there is a SZ fitting
	FI0[i]=Fi_initial_helix_referenceframe[i];

	// flag of a good fit in SZ space; later they will be set at true if fit in SZ is successful;
	GoodSkewFit[i]=false;
	fnSttSkewHitsinTrack[i]=0;

 }  // end for(i=0; i<nSttTrackCand;i++)


//-------------- stampa
 if(istampa>=2){
	cout<<"from PndTrkTracking, before Ordering, after primo loop sui paralleli.\n";
	cout<<"\tevt. n. "<<IVOLTE<<", n. di cicli fatti nel primo loop "<<iconta<<endl;
	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,-1,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }

//-------------- fine stampa


// now the ordering the parallel and skew hits.


 for(i=0; i<nSttTrackCand;i++){

	// the BigList array must be loaded also when there are no Skew hits.
	OrderingSttSkewandSttParallel(
		fOx[i],
		fOy[i],
		fR[i],
		fnSttSkewHitsinTrack[i],  // input
		&fListSttSkewHitsinTrack[i][0],// input, but this gets ordered
		&Sfinal[i][0], // input from Skew Straws
		Charge[i], // input
		fnSttParHitsinTrack[i], // input, # parallel Hits in the current track
		&fListSttParHitsinTrack[i][0],	// this was already ordered
		&U[i][0], // U conformal parallel hits; input, this was already ordered
		&V[i][0], // V conformal parallel hits; input, this was already ordered
		&BigList[i][0] // this is the final ordered Parallel+Skew list; already
				// in NATIVE hit number.
				);


 }   //  end of     for(i=0; i<nSttTrackCand;i++)

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


 nTotalCandidates = nSttTrackCand;  // nSttTrackCand is already <= MAXTRACKSPEREVENT.


//----   find the angular range (in Fi) allowed for the STT hits, with the present fOx,fOy and fR
//	of the track candidates, and for the Mvd hits (FI0 and Fi_low_limit ).

 for(  ncand= 0; ncand< nSttTrackCand; ncand++){

//------------------------------------------
fnMvdPixelHitsinTrack[ncand] =0;
fnMvdStripHitsinTrack[ncand] =0;
//--------------------------------------------

	keepit[ncand]=true;
	FindingParallelTrackAngularRange(
		fOx[ncand],
		fOy[ncand],
		fR[ncand],
		Charge[ncand],
		&Fi_low_limit[ncand],	// Fi (in XY Helix frame) lower limit using
					// the Stt detector minimum/maximum radius
					// Fi_low_limit is ALWAYS between 0. and 2PI
		&Fi_up_limit[ncand],	// Fi (in XY Helix frame) upper limit using
					// the Stt detector maximum/minimum radius
					// Fi_up_limit is ALWAYS > Fi_low_limit and
					// possibly > 2PI.
		&statusflag[ncand],//it is a vector; =0, all well; =1, track contained completely between RMin
			   // and RMax; = -1 track contained within RMin; =-2 track outside RMax.
		RSTRAWDETECTORMIN,
		RSTRAWDETECTORMAX
		);

	if( statusflag[ncand] == -1) { Fi_low_limit[ncand] = -99999.; }
	else if (statusflag[ncand] == -2){ keepit[ncand] = false; };

 }	// end of for(  ncand= 0; ncand< nSttTrackCand; ncand++)

//---------------------   here call to the function that matches Mvd hits with Stt tracks
   delta=0.5; //  parameter of proximity for associating Mvd hits to Stt tracks
   highqualitycut=0.2; //  parameter of proximity for associating Mvd hits to Stt tracks


   //  This method matches the Mvd hits to the found tracks.

   MatchMvdHitsToSttTracks(
			keepit,// input and output.
			delta,
			highqualitycut,
			nSttTrackCand,
			FI0,
			Fi_low_limit,
			Charge,
			fnMvdPixelHitsinTrack, // output
			fListMvdPixelHitsinTrack, // output
			fnMvdStripHitsinTrack, // output
			fListMvdStripHitsinTrack // output
			);

//-------------- stampa
 if(istampa>=2){
	cout<<"dopo MatchMvdHitsToSttTracks[1646], evt. "<<IVOLTE<<endl;
	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,-1,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa


//-----------------  end of section with match Mvd hits with Stt hits

	for(  ncand= 0; ncand< nSttTrackCand; ncand++){

		if(!keepit[ncand]) continue;

		if(fnMvdPixelHitsinTrack[ncand]+
			fnMvdStripHitsinTrack[ncand]==0){

		   Mvdhits[ncand]=false;
		   fALFA[ncand]=-2.*fOx[ncand];
		   fBETA[ncand]=-2.*fOy[ncand];
		   fGAMMA[ncand]= fOx[ncand]*fOx[ncand]+fOy[ncand]*fOy[ncand]-fR[ncand]*fR[ncand];
		   continue;
		} // end of  if(fnMvdPixelHitsinTrack[ncand]+


		if( fnMvdPixelHitsinTrack[ncand] > MAXMVDPIXELHITSINTRACK){
			cout<<"from PndTrkTracking, fnMvdPixelHitsinTrack["<<ncand
			<<"] is "<<fnMvdPixelHitsinTrack[ncand]
			<<"and it is > MAXMVDPIXELHITSINTRACK ("
			<<MAXMVDPIXELHITSINTRACK
			<<"); setting fnMvdPixelHitsinTrack[ncand] to "<<
			MAXMVDPIXELHITSINTRACK<<endl;
			fnMvdPixelHitsinTrack[ncand]=MAXMVDPIXELHITSINTRACK;
		}

		if( fnMvdStripHitsinTrack[ncand] > MAXMVDSTRIPHITSINTRACK){
			cout<<"from PndTrkTracking, fnMvdStripHitsinTrack["<<ncand
			<<"] is "<<fnMvdStripHitsinTrack[ncand]
			<<"and it is > MAXMVDSTRIPHITSINTRACK ("
			<<MAXMVDSTRIPHITSINTRACK
			<<"); setting fnMvdStripHitsinTrack[ncand] to "<<
			MAXMVDSTRIPHITSINTRACK<<endl;
			fnMvdStripHitsinTrack[ncand]=MAXMVDSTRIPHITSINTRACK;
		}


		Mvdhits[ncand] = true;



//----------------------------------------------------------
//     ordering all the hits belonging to the candidate track, by increasing fR;
//     forming the new track with Mvd+Stt hits


		// here fnTrackCandHit[ncand] is the sum of all the hits (Mvd+Stt) belonging to the track;
		fnTrackCandHit[ncand] =fnSttParHitsinTrack[ncand]+
					fnSttSkewHitsinTrack[ncand]+
					fnMvdPixelHitsinTrack[ncand]+
					fnMvdStripHitsinTrack[ncand];
	// arrays used to store temporarily the info of Mvd hits to be ordered.
		// adding the Mvd hits (Pixel and Strips)
			for(i=0; i< fnMvdPixelHitsinTrack[ncand]; i++){
			  XY[i][0] = fXMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ];
			  XY[i][1] = fYMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ];
			  ListHits[i] = fListMvdPixelHitsinTrack[ncand][i];
			}
			for(i=0; i< fnMvdStripHitsinTrack[ncand]; i++){
			  XY[i+fnMvdPixelHitsinTrack[ncand]][0] =
			    fXMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ];
			  XY[i+fnMvdPixelHitsinTrack[ncand]][1] =
			    fYMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ];
			  // to distinguish between Pixels and Strips, add a number
			  // to the original Strip hit number.
			  ListHits[i+fnMvdPixelHitsinTrack[ncand]] =
				fListMvdStripHitsinTrack[ncand][i]+
				(MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*10 ;
			}

			//  ordering the Mvd Hits
			OrderingUsingConformal(
			   fOx[ncand],
			   fOy[ncand],
			   &Trajectory_Start[ncand][0],
			   fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand],
			   XY, // XY[*][0] = X position, XY[*][0] = Y position.
			   Charge[ncand],  // input
			   tListHits  // output
						);

			//  constructing the ordered new Track  Candidate now and loading
			//  the (now ordered) X and Y lists for the hits in this track.
			for(i=0; i< fnMvdPixelHitsinTrack[ncand]+
					fnMvdStripHitsinTrack[ncand]; i++){
				if(ListHits[i]<(MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*10){//Pixel.
					fListTrackCandHit[ncand][i] = ListHits[i];
					fListTrackCandHitType[ncand][i] = 0;
				} else { // Strip hits.
					fListTrackCandHit[ncand][i] = ListHits[i]-
						(MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*10;
					fListTrackCandHitType[ncand][i] = 1;
				}
			}  // end of  for(i=0; i< fnMvdPixelHitsinTrack[ncand]+



		for(i=0; i<fnSttParHitsinTrack[ncand]+
			fnSttSkewHitsinTrack[ncand];i++){
			fListTrackCandHit[ncand][i+
				fnMvdPixelHitsinTrack[ncand]+
				fnMvdStripHitsinTrack[ncand]
					] = BigList[ncand][i];

			if( info[BigList[ncand][i]][5] < 2.){
				// parallel Stt hit.
				fListTrackCandHitType[ncand][i+
				fnMvdPixelHitsinTrack[ncand]+
				fnMvdStripHitsinTrack[ncand]
					] = 2;
			} else {
				fListTrackCandHitType[ncand][i+
				fnMvdPixelHitsinTrack[ncand]+
				fnMvdStripHitsinTrack[ncand]
					] = 3;
			}
		}


	} //   end of  for(ncand=0; ncand< nTotalCandidates; ncand++)


//-------------------- end of ordering

//-------------- stampa
 if(istampa>=2){
	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,-1,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa


//---------------------   refit the Helix in XY plane using Stt + Mvd associated hits


	Short_t	iexcl;
	Double_t
		d,
		diff,
		rotationangle;




 for(ncand=0; ncand< nTotalCandidates; ncand++){
	if(!keepit[ncand]) continue;

	if( !Mvdhits[ncand]) continue;


	//   Py/Px = m in    v = m*u + q formula
	rotationangle= atan2( Charge[ncand]*fOx[ncand],-Charge[ncand]*fOy[ncand]);	// Py/Px = m in
	if( rotationangle<0. ) rotationangle += 2.*PI;


//---------- translation of the reference system in the best Mvd hit position

	for(i=0, diff = 1.e20;i<fnTrackCandHit[ncand];i++){
// fListTrackCandHitType[ncand][i]=0 --> Pixel; = 1 --> Strip; = 2 Stt parallel; 3 = Stt Skew.
// There are no -1 type hits at this point.
		if(!(fListTrackCandHitType[ncand][i]==1
				||
				fListTrackCandHitType[ncand][i]==0) ) continue;
		if( fListTrackCandHitType[ncand][i] == 0 ){
				d = fabs( sqrt(
				(fOx[ncand]-fXMvdPixel[fListTrackCandHit[ncand][i]])
				*(fOx[ncand]-fXMvdPixel[fListTrackCandHit[ncand][i]])
				+(fOy[ncand]-fYMvdPixel[fListTrackCandHit[ncand][i]])
				*(fOy[ncand]-fYMvdPixel[fListTrackCandHit[ncand][i]])
					) - fR[ncand]);
			if( d  <  diff ){
				diff = d;
				iexcl=i;
				trajectory_vertex[0]=
					fXMvdPixel[fListTrackCandHit[ncand][i]];
				trajectory_vertex[1]=
					fYMvdPixel[fListTrackCandHit[ncand][i]];
			}
		} else {
			d = fabs( sqrt(
				(fOx[ncand]-fXMvdStrip[fListTrackCandHit[ncand][i]])
			       *(fOx[ncand]-fXMvdStrip[fListTrackCandHit[ncand][i]])
			       +(fOy[ncand]-fYMvdStrip[fListTrackCandHit[ncand][i]])
			       *(fOy[ncand]-fYMvdStrip[fListTrackCandHit[ncand][i]])
					) - fR[ncand]);
			if( d  <  diff ){
				diff = d;
				iexcl=i;
				trajectory_vertex[0]=
						fXMvdStrip[fListTrackCandHit[ncand][i]];
				trajectory_vertex[1]=
						fYMvdStrip[fListTrackCandHit[ncand][i]];
			}
		}

	} // end of  for(i=0, diff = 1.e20;i<.....

	if( diff > 0.5 ) {
		trajectory_vertex[0]=trajectory_vertex[1]=0.;
		iexcl=-1;
	}

//---------------- inizio stampe.
if(istampa>=2){
	if(iexcl!= -1){
		cout<<"from PndSttMvdTracking, evt. "<<IVOLTE<<", cand n. "<<ncand
		<<", Prima di refit, Pivoting hit :\n";
		if(fListTrackCandHitType[ncand][iexcl]==0){
			cout<<"\t Pixel hit n. "<<fListTrackCandHit[ncand][iexcl]<<endl;
		} else {
			cout<<"\t Strip hit n. "<<fListTrackCandHit[ncand][iexcl]<<endl;
		}
	} else {
		cout<<"from PndSttMvdTracking, evt. "<<IVOLTE
			<<", Prima di refit, no Pivoting hit.\n";
	}
	cout<<"\trotation angle = "<<rotationangle<<endl;
}
//-------------- fine di stampe.

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

//	trajectory_vertex[0]=0.;
//	trajectory_vertex[1]=0.;
//	iexcl = -1;
//	rotationangle=0.;


	RefitMvdStt(
			fnTrackCandHit[ncand], //this is  input
			&fListTrackCandHit[ncand][0],//this is both input and output
			&fListTrackCandHitType[ncand][0],//this is both input and output
			info,
			rotationangle,
			trajectory_vertex,
			iexcl,

			&fALFA[ncand], // output of the fit
			&fBETA[ncand], // output of the fit
			&fGAMMA[ncand],// output of the fit
			&status[ncand]    // fit status; true = successful
			);

	if(status[ncand]){
			fOx[ncand] = -fALFA[ncand]/2.;
			fOy[ncand] = -fBETA[ncand]/2.;
			fR[ncand] = fOx[ncand]*fOx[ncand]+fOy[ncand]*fOy[ncand]-fGAMMA[ncand];
			if( fR[ncand] > 0. ) fR[ncand]=sqrt(fR[ncand]) ;
			else  fR[ncand]=0.;
			FI0[ncand] = atan2(Trajectory_Start[1][ncand]-fOy[ncand], Trajectory_Start[0][ncand]-fOx[ncand]);
			if( FI0[ncand] < 0. )  FI0[ncand]+= 2.*PI;
	}

	if(fR[ncand]<RSTRAWDETECTORMIN/2.) {
		keepit[ncand]=false;
		continue;
	}
//--------------

//  find out if the SciTil hit associated with this track is still acceptable after the last
//  XY refit.


	for(i=0, iaccept=0;i<fnSciTilHitsinTrack[ncand];i++){

		intersect=GeomCalculator.IntersectionSciTil_Circle(
			DIMENSIONSCITIL,
			fposizSciTil[fListSciTilHitsinTrack[ncand][i]][0],
			fposizSciTil[fListSciTilHitsinTrack[ncand][i]][1],
			fOx[ncand], // center of circle.
			fOy[ncand],
			fR[ncand], // Radius of circle.
			&Nint,
			XintersectionList,
			YintersectionList
						);

		// reject case with no intersection of the SciTil with the circle trajectory.
		if(intersect){
		// calculate S on the lateral face of the Helix.
			if ( Nint==2){// (it should be a very rare case).
				// do an average of the two positions; put it in XintersectionList[0]
				// and in  YintersectionList[0] so that it is transparent to
				// the subsequent code.
				YintersectionList[0]=
					(YintersectionList[0]+YintersectionList[1])/2.;
				XintersectionList[0]=
					(XintersectionList[0]+XintersectionList[1])/2.;
			} // end of  if ( Nint==2)
			fListSciTilHitsinTrack[ncand][iaccept]= 
					fListSciTilHitsinTrack[ncand][i];
			fS_SciTilHitsinTrack[ncand][0]=atan2(YintersectionList[0]-fOy[ncand],
				XintersectionList[0]-fOx[ncand]);
			if ( fS_SciTilHitsinTrack[ncand][0]<0.)
				fS_SciTilHitsinTrack[ncand][0] += 2.*PI;
			iaccept++;


		}   //  end of  if(intersect)

	} // end of  for(i=0, iaccept=0;i<fnSciTilHitsinTrack[ncand];i++)



	fnSciTilHitsinTrack[ncand]=iaccept;



// finding again the Charge of the track (the last fit may have changed the concavity of the track).

// first load again the (now ordered) X and Y position of the hits.

	for(i=0;i<fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand];i++){
		// these are not ordered, but for charge finding ordering is not
		// necessary.
		X[i]=XY[i][0];
		Y[i]=XY[i][1];
	}

	j = fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand];
	for(i=0;i<fnSttParHitsinTrack[ncand];i++){
		X[j]=info[ fListSttParHitsinTrack[ncand][i] ][0];
		Y[j]=info[ fListSttParHitsinTrack[ncand][i] ][1];
		j++;
	}
	for(i=0;i<fnSciTilHitsinTrack[ncand];i++){
		X[j]=fposizSciTil[ fListSciTilHitsinTrack[ncand][i] ][0];
		Y[j]=fposizSciTil[ fListSciTilHitsinTrack[ncand][i] ][1];
		j++;
	}

	FindCharge(
		fOx[ncand],
		fOy[ncand],
		fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]
		   +fnSttParHitsinTrack[ncand]+fnSciTilHitsinTrack[ncand],
		X,
		Y,
		&Charge[ncand]
		);

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


	//   finding the FI angular range (in the laboratory frame) spanned by this parallel track
	//   taking into account the Rmax and Rmin of the straw detector. The track NOT NECESSARILY
	//   has to originate at (0,0).

	FindingParallelTrackAngularRange(
		fOx[ncand],
		fOy[ncand],
		fR[ncand],
		Charge[ncand],  // this is an input.
		&Fi_low_limit[ncand],	// Fi (in XY Helix frame) lower limit using
					// the Stt detector minimum/maximum radius
					// Fi_low_limit is ALWAYS between 0. and 2PI
		&Fi_up_limit[ncand],	// Fi (in XY Helix frame) upper limit using
					// the Stt detector maximum/minimum radius
					// Fi_up_limit is ALWAYS > Fi_low_limit and
					// possibly > 2PI.
					// Fi_low_limit and Fi_up_limit are undefined
					// when statusflag is negative.
		&statusflag[ncand],//it is a vector; =0, all well; =1, track contained completely between RMin
			// and RMax; = -1 track contained within RMin; =-2 track outside RMax.
		RSTRAWDETECTORMIN,
		RSTRAWDETECTORMAX
						);
	// Fi_low_limit set  to  -99999. when it is contained in Mvd
	// region completely.
	if( statusflag[ncand] == -1 ) Fi_low_limit[ncand] = -99999.;
	if( statusflag[ncand] == -2) keepit[ncand] = false;


	}	// end of for(ncand=0; ncand< nTotalCandidates; ncand++)

//--------------------- end of  refit the Helix in XY plane using Stt + Mvd associated hits


//-------------- stampa
 if(istampa>=2){
	cout<<"\tstampa dopo il Refit.\n";

	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,-1,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa


  if(fnMvdPixelHit+fnMvdStripHit>0){
//---------------------   here call to the function that matches Mvd hits with the new
//  circular trajectory in XY found  for the second time, after first refit

   delta=0.5; //  parameter of proximity for associating Mvd hits to Stt tracks
//   highqualitycut=0.3; //  parameter of proximity for associating Mvd hits to Stt tracks
   highqualitycut=0.5; //  parameter of proximity for associating Mvd hits to Stt tracks

   MatchMvdHitsToSttTracksagain(
			keepit,
			Mvdhits,
			delta,
			highqualitycut,
			nSttTrackCand,
			FI0,
			Fi_low_limit, // because here we deal with hits in Mvd region
			Charge,
			fnMvdPixelHitsinTrack, // input and output
			fListMvdPixelHitsinTrack, // input and output
			fnMvdStripHitsinTrack, // input and output
			fListMvdStripHitsinTrack // input and output
			);
  } // end of  if(fnMvdPixelHit+fnMvdStripHit>0)

//---------------------
// check if the track candidate has still at least 2 axial hits;
  for(ncand=0; ncand< nTotalCandidates; ncand++)
  {
	if( fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+fnSttParHitsinTrack[ncand]<2) keepit[ncand]=false;
  }  // end of for(ncand=0; ncand< nTotalCandidates; ncand++)
//---------------------   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


//-------------- stampa
 if(istampa>=2){
cout<<"\tstampa dopo il Match Again.\n";
fPrint.stampetta(IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],
&fListMvdStripHitsinTrack[0][0],&fListSttParHitsinTrack[0][0],
&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,
fnSttSkewHitsinTrack,fnSciTilHitsinTrack,nSttTrackCand,
-1, // print all candidates;
MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa



//-----------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------
//----------- from here on, the XY parameters of the Helix don't change. ------------------
//-----------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------

// use the risult just obtained from the fit in XY to do the association of the Skew Straw hits


  for(ncand=0; ncand< nTotalCandidates; ncand++)
  {
	if(!keepit[ncand]) continue;

	if( statusflag[ncand] == -1 ) continue; // this is when the XY circle is contained in the
						// the Mvd region completely; skip the association of
						// the Skews.

	// the arrays  Temporary....   are related to the current candidate; they loose meaning
	// as soon as the loop over the candidate tracks finishes;

	fnSttSkewHitsinTrack[ncand]= AssociateSkewHitsToXYTrack(
		fInclusionListStt, // hit is excluded only if it multiple hit.
		nSttSkewHit,
		fListSttSkewHits,
		fOx[ncand],   //  input : X of center of XY plane circle
		fOy[ncand],   //  input : Y of center of XY plane circle
		fR[ncand],   //  input : Radius of XY plane circle
		info,
		WDX,
		WDY,
		WDZ,
		Fi_low_limit[ncand],	// Fi (in Helix XY frame) lower limit using the Stt detector minimum/maximum radius
		Fi_up_limit[ncand],	// Fi (in Helix XY frame) upper limit using the Stt detector maximum/minimum radius
		Charge[ncand],

		TemporarySkewList, // output : TemporarySkewList[*][0] =  list of selected skew hits
				   // numbers(in skew numbering);
				   // TemporarySkewList[*][1] =  solution number (0 or 1) as per
				   // the calculateintersections  method;
	// the info on TemporaryS, TemporaryZ, TemporaryZDrift, TemporaryZErrorafterTilt  as follows :
	// if   i  is the (original) Skew Hit number, and ii (0 or 1) is
	// the solution, then the infos are stored in the   i + ii*MAXSTTHITS location;
		TemporaryS,       //  output,  S coordinate of selected Skew hit
		TemporaryZ,       //  output,  Z coordinate of selected Skew hit (center wire)
		TemporaryZDrift,   //  output,  drift distance IN Z DIRECTION only, of selected Skew hit
		TemporaryZErrorafterTilt   //  output,  Radius taking into account the tilt, IN Z DIRECTION only, of selected Skew hit
			);



	// limit the total # Stt hits to MAXSTTHITSINTRACK
	if( fnSttSkewHitsinTrack[ncand]+fnSttParHitsinTrack[ncand] > MAXSTTHITSINTRACK ) {
	  if(MAXSTTHITSINTRACK > fnSttParHitsinTrack[ncand])
		fnSttSkewHitsinTrack[ncand]=MAXSTTHITSINTRACK-fnSttParHitsinTrack[ncand];
	  else fnSttSkewHitsinTrack[ncand]=0;
	}
	for(j=0;j<fnSttSkewHitsinTrack[ncand];j++)
	{
		fListSttSkewHitsinTrack[ncand][j]=TemporarySkewList[j][0];
		fListSttSkewHitsinTrackSolution[ncand][j]=TemporarySkewList[j][1];
	}

//-------------------------------------------  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//-------------- stampa
 if(istampa>=2){
cout<<"\tstampa dopo AssociateSkewHitsToXYTrack;\n";fPrint.stampetta(IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],
&fListMvdStripHitsinTrack[0][0],&fListSttParHitsinTrack[0][0],
&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,
fnSttSkewHitsinTrack,fnSciTilHitsinTrack,nSttTrackCand,
ncand, // print the candidate  ncand;
MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa


// calculate the number of hits to use in the SZ fit later;
	//  nXYZhits = n. of Mvd hits + SciTil hits. However, if there are 2 SciTil
	//  hits in this track (namely two adjacent SciTil tiles have a hit
	//  caused PRESUMABLY by the same track) then count them AS ONE because below
	//  the average of their postions is considered !

	if( fnSciTilHitsinTrack[ncand] == 2) {
	   nXYZhits = fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+ 1;
	}else{   // in this case fnSciTilHitsinTrack[ncand] is 0 or 1;
	   nXYZhits = fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+
			fnSciTilHitsinTrack[ncand];
	}

	// calculate if there is the need of using some skew hits in the subsequent SZ fit;
	// put in  nhitsinfit  the number of hits used in the subsequent  SZ  fit.


	// the following is valid only for GLPK fits;
	if(YesGLPKfitSZ)
	{	// flag for the GLPK choice of fit;
	  if( nXYZhits <=2){
		fnSttSkewHitsinTrack[ncand]<5 ?
			nhitsinfit = nXYZhits + fnSttSkewHitsinTrack[ncand] :
			nhitsinfit = nXYZhits + 5 ; // 1 0 2 XYZ hit + 5 Skew hits.
	  } else {
		nhitsinfit= nXYZhits;
	  }
	} else {  // other fit choice;

	  // the following is valid for non-GLPK fits;
	  nhitsinfit = nXYZhits + fnSttSkewHitsinTrack[ncand];
	}  // end of if(YesGLPKfitSZ)

	// the following is a protection against declaration of 0 dimension array;
	int dime ;
	if(nhitsinfit>0) dime = nhitsinfit ; else dime=1;


//---------------------   here calculate the S and Z values of Mvd Pixels, Mvd Strips,
//	 Stt Skew hits and SciTil hits (if present).

//  the difference between S and Sbis, ZED and ZEDbis, DriftRadius and DriftRadiusbis,
//  ErrorDriftRadius and ErrorDriftRadiusbis, is that S, ZED etc. contain the list of
//  Pixel+Strips+SciTil + other Skew Stt hits in case Pixel+Strips+SciTil are <= 2; instead
//  Sbis, ZEDbis etc. contain Pixel+Strips+all Skew Stt hits.


	// load the quantities needed for the SZ fit;

	LoadSZetc_forSZfit(
		ncand,	// input
		nhitsinfit,
		TemporaryS,		// input
		TemporaryZ,		// input
		TemporaryZDrift,	// input
		TemporaryZErrorafterTilt,	// input
		YesGLPKfitSZ,		// input

		ErrorDriftRadius,	 // output
		ErrorDriftRadiusbis,	 // output
		DriftRadius,		 // output
		DriftRadiusbis,	 // output
		S,			 // output
		Sbis,		 // output
		ZED,			 // output
		ZEDbis		 // output
	);


// ---------------  fit in SZ  with Mvd + SciTil only



//  finding if there are discontinuity at 0 for fi value of the Mvd Hit.
//  In case of discontinuity at 0, add 2*PI to fi of those hits with fi in the 1st quadrant.
//  This is necessary because the discontinuities would make the fit
//  in the SZ plane fail.
//  In this discontinuity fixing, the value FI0 of the vertex (0,0) is also included.
//  If there is discontinuity fixing, the values of S[i] might be modified (+2.*PI).
		FixDiscontinuitiesFiangleinSZplane(
			nhitsinfit,
			S,	// S can be modified by +-2*PI if necessary.
			&FI0[ncand],	// this remains unchanged.
			Charge[ncand]	// this remains unchanged.
					);

//-------------- stampa
 if(istampa>=2){
	cout<<"\tstampa prima di FitSZspace, [non in Mvd track section] .\n";
	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,ncand,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);

 }
//-------------- fine stampa
//---------------------   here do the fit again in the SZ space if there are Mvd hits.
//			  For this, reordering of the  Mvd hits is not necessary.


	if(nhitsinfit>0){
		resultFitSZagain[ncand] = fit.FitSZspace(
				nhitsinfit,	// n. hits to be fitted
				tS,
				tZED,
				tDriftRadius,
				tErrorDriftRadius,
				FI0[ncand],
				MAXSKEWHITSINFIT,// maximum number of STT Skew hits in fit;
				&emme,
				IVOLTE*100+ncand // number of the accumulation plot.
						);


//-------------- stampa
 if(istampa>=2){
	cout<<"\tstampa dopo FitSZspace, [not in Mvd track section], result (1 va bene) = "
	<<resultFitSZagain[ncand]<<endl;
	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,ncand,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
	}
//---------------------------------------------------------------------------
		if( resultFitSZagain[ncand]==1){
			KAPPA[ncand] = emme;
			GoodSkewFit[ncand] = true;
			if( ncand<= nSttTrackCand ) SttSZfit[ncand]=true;
		} else {
			keepit[ncand]=false;
			GoodSkewFit[ncand] = false;
		}


	}  else {  // continuation of  if(nhitsinfit>0)
		keepit[ncand]=false;
		GoodSkewFit[ncand] = false;
	}   // end of  if(nhitsinfit>0)



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

//	use the risult just obtained from the fit in SZ to reject the spurious Skew Straw hits
//	and the Mvd spurious hits; also in this function there is the calculation of the Z position
//	of the SKEW hits and the MVD hits, for a given track candidate (ie for a given Helix
//	circle in the XY plane)

	Short_t MaxTurns;
	Double_t Turns;


	  if(keepit[ncand]){
	    if(fR[ncand] < RSTRAWDETECTORMAX/2.){
		if(-Charge[ncand]*KAPPA[ncand]>0.){	// this means Pz>0.
		  Turns= 0.5*fabs((ZCENTER_STRAIGHT+SEMILENGTH_STRAIGHT)
					*KAPPA[ncand])/PI;
		  if( fabs(Turns)<10.)  MaxTurns=(Short_t) Turns ;
		  else  MaxTurns=10;
		} else {
		  Turns= 0.5*fabs((ZCENTER_STRAIGHT-SEMILENGTH_STRAIGHT)
					*KAPPA[ncand])/PI;
		  if( fabs(Turns)<10.)  MaxTurns=(Short_t) Turns ;
		  else  MaxTurns=10;
		}
	    } else {
		MaxTurns=0;
	    }

	    EliminateSpuriousSZ(
		MaxTurns,
		&fnMvdPixelHitsinTrack[ncand],	// input and output
		&fListMvdPixelHitsinTrack[ncand][0],// input and output
		&fnMvdStripHitsinTrack[ncand],	// input and output
		&fListMvdStripHitsinTrack[ncand][0],// input and output
		&fnSttSkewHitsinTrack[ncand],	// input and output
		&fListSttSkewHitsinTrack[ncand][0],// input and output
		Sbis,	// input, position of the central wire on the Helix cylinder;
		ZEDbis,	// input, position of the central wire on the Helix cylinder.
		DriftRadiusbis,	// input
		ErrorDriftRadiusbis,	// input
		&SchosenPixel[ncand][0], // this value from now on
		&SchosenStrip[ncand][0], // can also be > 2PI or < 2PI when
		&SchosenSkew[ncand][0],  // the particle makes more than 1 turn.
		&ZchosenPixel[ncand][0],
		&ZchosenStrip[ncand][0],
		&ZchosenSkew[ncand][0],
		ErrorchosenPixel,
		ErrorchosenStrip,
		ErrorchosenSkew,
		KAPPA[ncand],
		FI0[ncand],
		fR[ncand]
		    );

	  }  // end of  if(keepit[ncand])

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


//-------------- stampa
 if(istampa>=2){
	cout<<"\tstampa dopo EliminateSpuriousSZ, [non in Mvd track section] .\n";

	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,ncand,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa

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

// redo fit after spurious cleaning;
		// recalculate the input values with the new list of hits
		// belonging to this track cand;

/*
		resultFitSZagain[ncand] = fit.FitSZspace(
				nhitsinfit,	// n. hits to be fitted
				tS,
				tZED,
				tDriftRadius,
				tErrorDriftRadius,
				FI0[ncand],
				MAXSKEWHITSINFIT,// maximum number of STT Skew hits in fit;
				&emme,
				0 // number of the accumulation plot;no meaning
						// here for this type of fit;
						);
		if( resultFitSZagain[ncand]==1){
			KAPPA[ncand] = emme;
			GoodSkewFit[ncand] = true;
			if( ncand<= nSttTrackCand ) SttSZfit[ncand]=true;
		} else {
			keepit[ncand]=false;
			GoodSkewFit[ncand] = false;
		}


//    redo the spurious cleaning after the second iteration SZ fit;

	  if(keepit[ncand]){
	    if(fR[ncand] < RSTRAWDETECTORMAX/2.){
		if(-Charge[ncand]*KAPPA[ncand]>0.){	// this means Pz>0.
		  Turns= 0.5*fabs((ZCENTER_STRAIGHT+SEMILENGTH_STRAIGHT)
					*KAPPA[ncand])/PI;
		  if( fabs(Turns)<10.)  MaxTurns=(Short_t) Turns ;
		  else  MaxTurns=10;
		} else {
		  Turns= 0.5*fabs((ZCENTER_STRAIGHT-SEMILENGTH_STRAIGHT)
					*KAPPA[ncand])/PI;
		  if( fabs(Turns)<10.)  MaxTurns=(Short_t) Turns ;
		  else  MaxTurns=10;
		}
	    } else {
		MaxTurns=0;
	    }

	    EliminateSpuriousSZ(
		MaxTurns,
		&fnMvdPixelHitsinTrack[ncand],	// input and output
		&fListMvdPixelHitsinTrack[ncand][0],// input and output
		&fnMvdStripHitsinTrack[ncand],	// input and output
		&fListMvdStripHitsinTrack[ncand][0],// input and output
		&fnSttSkewHitsinTrack[ncand],	// input and output
		&fListSttSkewHitsinTrack[ncand][0],// input and output
		Sbis,	// input, position of the central wire on the Helix cylinder;
		ZEDbis,	// input, position of the central wire on the Helix cylinder.
		DriftRadiusbis,	// input
		ErrorDriftRadiusbis,	// input
		&SchosenPixel[ncand][0], // this value from now on
		&SchosenStrip[ncand][0], // can also be > 2PI or < 2PI when
		&SchosenSkew[ncand][0],  // the particle makes more than 1 turn.
		&ZchosenPixel[ncand][0],
		&ZchosenStrip[ncand][0],
		&ZchosenSkew[ncand][0],
		ErrorchosenPixel,
		ErrorchosenStrip,
		ErrorchosenSkew,
		KAPPA[ncand],
		FI0[ncand],
		fR[ncand]
		    );

	  }  // end of  if(keepit[ncand])

*/
//------------------------

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


//	First cleanup based on the absence of Mvd hits


	if(fYesCleanMvd){

		// reject the candidate if it is NOT contained in the pipe and
		// therefore it should have at least 1 Mvd hit but it has none.
		if( (!GeomCalculator.IsInTargetPipe(
			fOx[ncand],
			fOy[ncand],
			fR[ncand],
			FI0[ncand],
			KAPPA[ncand],
			Charge[ncand],
			VERTICALGAP/2.) )
				 &&
			fnMvdStripHitsinTrack[ncand]+fnMvdPixelHitsinTrack[ncand]==0)
		{
			keepit[ncand]=false;
		}
	}  // end of  (fYesCleanMvd)

  }	//  end of for(ncand=0; ncand< nTotalCandidates; ncand++)


 // In case some candidate track were not processed in the previous loop, there might
 // be still GoodSkewFit[ncand] = false; in that case the candidate track is rejected
 // because it has no information on KAPPA;
  for(ncand=0; ncand< nTotalCandidates; ncand++){
	if(!GoodSkewFit[ncand]) keepit[ncand]=false;
  }

//--------

//-------------- stampa
 if(istampa>=2){
	cout<<"\tstampa dopo il Cleanup piccolo\n";
	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,-1,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa

//---  redo association of parallel Stt  straw  hits to this track, after better refit.
  if(fnMvdPixelHit+fnMvdStripHit>0){


	CollectParSttHitsagain(
			keepit,
			Mvdhits,
			info,
			nSttParHit,
			0,	// starting candidate n. (included).
			nTotalCandidates, // ending candidate n. (excluded).
			KAPPA,
			FI0,
			Fi_low_limit,
			Fi_up_limit,
			fnSttParHitsinTrack, // input and output
			fListSttParHitsinTrack // input and output
			);
//-------------------------------
  for(ncand=0; ncand< nTotalCandidates; ncand++){
	// check if the track candidate has still at least 2 axial hits;
	if( fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+fnSttParHitsinTrack[ncand]<2){
		keepit[ncand]=false;
		continue;
	}

	// limit the total # Stt hits to MAXSTTHITSINTRACK
	if( fnSttSkewHitsinTrack[ncand]+fnSttParHitsinTrack[ncand] > MAXSTTHITSINTRACK ) {
	  if(MAXSTTHITSINTRACK > fnSttSkewHitsinTrack[ncand])
		fnSttParHitsinTrack[ncand]=MAXSTTHITSINTRACK-fnSttSkewHitsinTrack[ncand];
	  else fnSttParHitsinTrack[ncand]=0;
	}
  }  // end of  for(ncand=0; ncand< nTotalCandidates; ncand++)
//---------------

//-------------- stampa
 if(istampa>=2){
	cout<<"\tstampa dopo CollectParSttHitsagain\n";
	fPrint.stampetta(
IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],&fListMvdStripHitsinTrack[0][0],
&fListSttParHitsinTrack[0][0],&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],
fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,fnSttSkewHitsinTrack,
fnSciTilHitsinTrack,nSttTrackCand,-1,MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,
MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa


  }  // end of if(fnMvdPixelHit+fnMvdStripHit>0
//--------------------- end of the section where the Mvd hits are added to the track candidates



//   the following is valid for all tracks;

//	ordering all the hits belonging to the candidate track, by increasing fR (large
//	trajectories)  or Conformal variables (better for small trajectories);
//	from candidate n. 0 to candidate n. nTotalCandidates-1; loading fListTrackCandHit.
//	the array ordered are :
//	fListTrackCandHit, fListTrackCandHitType, fListSttParHitsinTrack, fListSttSkewHitsinTrack
//	and also at the end the SciTil hit (if present) is added.


	Ordering_Loading_ListTrackCandHit(
		keepit,
		0,
		nTotalCandidates,
		info,
		Trajectory_Start,
		Charge,
		SchosenSkew
		);

	// adding at the end the SciTil hits (if present).


	for(ncand=0; ncand< nTotalCandidates; ncand++){
		if(!keepit[ncand]) continue;
		i=fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+
			fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand];
		for(j=0;j<fnSciTilHitsinTrack[ncand];j++) {
			fListTrackCandHit[ncand][i+j]=fListSciTilHitsinTrack[ncand][j];
			fListTrackCandHitType[ncand][i+j] = 1001;
		}
	}  // end of for(ncand=0; ncand< nTotalCandidates; ncand++)



//------------- cleanup section.


	Start[0]=0.;
	Start[1]=0.;
	Start[2]=0.;
	gap = (Double_t) (VERTICALGAP);

//-------------- stampa
 if(istampa>=2){
	cout<<"\tstampa prima del Cleanup grosso\n";
fPrint.stampetta(IVOLTE,tkeepit,&fListMvdPixelHitsinTrack[0][0],
&fListMvdStripHitsinTrack[0][0],&fListSttParHitsinTrack[0][0],
&fListSttSkewHitsinTrack[0][0],&fListSciTilHitsinTrack[0][0],fnMvdPixelHitsinTrack,fnMvdStripHitsinTrack,fnSttParHitsinTrack,
fnSttSkewHitsinTrack,fnSciTilHitsinTrack,nSttTrackCand,
-1, // print all candidates;
MAXMVDPIXELHITSINTRACK,MAXMVDSTRIPHITSINTRACK,MAXSCITILHITSINTRACK,MAXSTTHITSINTRACK,fR,fOx,fOy,FI0,KAPPA);
 }
//-------------- fine stampa


    for(ncand=0, nRemainingCandidates=0; ncand< nTotalCandidates; ncand++){
	if(!keepit[ncand]) continue;
	Short_t &nHitsPar = fnSttParHitsinTrack[ncand];
	Short_t &nHitsSkew = fnSttSkewHitsinTrack[ncand];
	Double_t auxS[nHitsSkew];

	for(i=0;i<nHitsSkew;i++){
		auxS[i] = SchosenSkew[ncand][fListSttSkewHitsinTrack[ncand][i]];
	}

	if(fYesClean){
		if(istampa>1) cout<<"PndTrkTracking, entra in TrackCleanup "<<
			"tracce normali, IVOLTE "<<IVOLTE<<" e track cand. "<<ncand<<endl;
		if ( !Cleaner.TrackCleanup(
		APOTEMAMAXINNERPARSTRAW,
		APOTEMAMAXSKEWSTRAW,
		APOTEMAMINOUTERPARSTRAW,
		APOTEMAMINSKEWSTRAW,
		auxS,
		Charge[ncand],
		FI0[ncand],
		gap,
		info,
		istampa,
		IVOLTE,
		KAPPA[ncand],
		&fListSttParHitsinTrack[ncand][0],
		&fListSttSkewHitsinTrack[ncand][0],
		MAXSTTHITS,
		nHitsPar,
		nHitsSkew,
		fOx[ncand],
		fOy[ncand],
		fR[ncand],
		RSTRAWDETECTORMAX,
		RSTRAWDETECTORMIN,
		SEMILENGTH_STRAIGHT,
		Start,
		STRAWRADIUS,
		ZCENTER_STRAIGHT
				) ) {
			keepit[ncand]=false;
			continue;
		} // end if

	}  // end of if(fYesClean)

	nRemainingCandidates++;

    }	//  end of for(ncand=0; ncand< nTotalCandidates; ncand++)

//------------- end of cleanup section.


//---------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------
//   this section collects more tracks starting from the remaining Mvd hits
//   not associated to any Stt track.

	if(fMvdAloneTracking && fnMvdPixelHit+fnMvdStripHit>0) {

//  load the UsedPixel and UsedStrip  vectors.
	bool	UsedPixel[MAXMVDPIXELHITS],
		UsedStrip[MAXMVDSTRIPHITS];

	Short_t List[MAXMVDPIXELHITS+MAXMVDSTRIPHITS],
		 ListType[MAXMVDPIXELHITS+MAXMVDSTRIPHITS];

	Double_t centrex,
		 centrey,
		 radius;

	for(i=0;i<fnMvdPixelHit; i++){
		UsedPixel[i]=false;
	}
	for(i=0;i<fnMvdStripHit; i++){
		UsedStrip[i]=false;
	}

	for(ncand=0; ncand< nTotalCandidates; ncand++){
//		if(!keepit[ncand]) continue;

		// here don't use the  keepit  information because
		// I don't want to risk to find the wrong tracks again.
		for(i=0; i< fnMvdPixelHitsinTrack[ncand]; i++){
			UsedPixel[ fListMvdPixelHitsinTrack[ncand][i] ] = true;
		}// end of for(i=0; i< fnMvdPixelHitsinTrack[ncand]; i++)

		for(i=0; i< fnMvdStripHitsinTrack[ncand]; i++){
			UsedStrip[ fListMvdStripHitsinTrack[ncand][i] ] = true;
		}// end of for(i=0; i< fnMvdStripHitsinTrack[ncand]; i++)

	}//  end of for(ncand=0; ncand< nTotalCandidates; ncand++)

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

	//  start from those Mvd CANDIDATES found by Riemann PR, whose MAJORITY of
	//  Pixel or Strip hits were not used previously in this code.

	for(i=0; i<fnMvdTrackCand; i++){
		if( nTotalCandidates >= MAXTRACKSPEREVENT -1 ) break; // protection for the
						// length of many arrays.
		nalone=0;
		npixelhitsintrack=0;
		nstriphitsintrack=0;
		for(j=0;j<fnHitMvdTrackCand[i];j++){

		  if(fListHitTypeMvdTrackCand[i][j] == 
			FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)){
		   if(UsedPixel[ fListHitMvdTrackCand[i][j] ])   continue;
		   if(npixelhitsintrack>MAXMVDPIXELHITSINTRACK) continue;
		   AloneX[nalone] = fXMvdPixel[ fListHitMvdTrackCand[i][j] ];
		   AloneY[nalone] = fYMvdPixel[ fListHitMvdTrackCand[i][j] ];
		   List[nalone] = fListHitMvdTrackCand[i][j];
		   ListType[nalone] = 0;
		   nalone++;
		   npixelhitsintrack++;
		  }else {// at this point this is a Strip hit; already made sure
			// earlier in the code that there is no third possibility.

//		   if (fListHitTypeMvdTrackCand[i][j] ==	// this condition in principle at this
//								// point is always satisfied.
//			FairRootManager::Instance()->GetBranchId(fMvdStripBranch)){


		   if(UsedStrip[ fListHitMvdTrackCand[i][j] ])   continue;
		   if(nstriphitsintrack>MAXMVDSTRIPHITSINTRACK) continue;
		   AloneX[nalone] = fXMvdStrip[ fListHitMvdTrackCand[i][j] ];
		   AloneY[nalone] = fYMvdStrip[ fListHitMvdTrackCand[i][j] ];
		   List[nalone] = fListHitMvdTrackCand[i][j];
		   ListType[nalone] = 1;
		   nalone++;
		   nstriphitsintrack++;
		  }
		}	// end of   for(j=0;j<fnHitMvdTrackCand[i];j++)

		if(nalone < 3 ) continue;

		// here find possibly new tracks

		if( nalone == 3 ){
			status[nTotalCandidates] = GeomCalculator.CalculateCircleThru3Points(
							AloneX[0],
							AloneY[0],
							AloneX[1],
							AloneY[1],
							AloneX[2],
							AloneY[2],
							&centrex,
							&centrey,
							&radius
							   );
			// the election on the radius is justified by the assumption that
			// the track originates from (0,0,0).
			if( status[nTotalCandidates] && radius > RSTRAWDETECTORMIN/2.){
				fOx[nTotalCandidates] = centrex;
				fOy[nTotalCandidates] = centrey;
				fR[nTotalCandidates] = radius;
				FI0[nTotalCandidates] =
				    atan2(-fOy[nTotalCandidates], -fOx[nTotalCandidates]);
				if( FI0[nTotalCandidates] < 0. )
					FI0[nTotalCandidates]+= 2.*PI;
				fnSttParHitsinTrack[nTotalCandidates]=0;
				fnSttSkewHitsinTrack[nTotalCandidates]=0;
				fnMvdPixelHitsinTrack[nTotalCandidates]=0;
				fnMvdStripHitsinTrack[nTotalCandidates]=0;
				for(j=0; j<nalone;j++){
				   if(ListType[j]==0){
					fListMvdPixelHitsinTrack[nTotalCandidates][
						fnMvdPixelHitsinTrack[nTotalCandidates]
										]=List[j];
						fnMvdPixelHitsinTrack[nTotalCandidates]++;
				   }else{
					fListMvdStripHitsinTrack[nTotalCandidates][
						fnMvdStripHitsinTrack[nTotalCandidates]
										]=List[j];
					fnMvdStripHitsinTrack[nTotalCandidates]++;
				   }
				}  //   end of  for(j=1; j<nalone;j++)

				fnTrackCandHit[nTotalCandidates]=nalone;
				keepit[nTotalCandidates]=true;
				Mvdhits[nTotalCandidates]=true;
				nTotalCandidates++;
				nRemainingCandidates++;
			} // end  of  if( status[nTotalCandidates] )


		} else { // continuation of  if( nalone == 3 )
			iexcl = 0;	// exclude first point in List from fit
					//  because it is taken as center of new
					//  reference frame in the fit.
			trajectory_vertex[0] = AloneX[0];
			trajectory_vertex[1] = AloneY[0];

			for(j=0; j<nalone;j++){
				fListTrackCandHit[nTotalCandidates][j]= List[j];
				fListTrackCandHitType[nTotalCandidates][j]= ListType[j];
			}
			// it is necessary to calculate rotationangle at this point [since this
			// is not an old stt candidate for which it is just enough to say
			// tan(rotationangle) = Py/Px ].  Use (almost) the same piece of code
			// existing in PndTrackFinderReal, modified by taking into account the
			// translation of axis in the Pivoting point.
			double rotationcos, rotationsin;
			for(j=1, rotationcos=0., rotationsin=0.; j<nalone; j++){
				// fi in the XY space is the same as fi in conformal space;
				double fi = atan2(AloneY[j]-AloneY[0],AloneX[j]-AloneX[0]);
				// use the cos and sin to avoid strange situations in the average when
				// the hits are aroung 0 degrees or 180 degrees;
				rotationcos += cos(fi);
				rotationsin += sin(fi);
			}
			rotationcos /= (nalone-1);
			rotationsin /= (nalone-1);
			rotationangle = atan2(rotationsin, rotationcos);
//----------------
if(istampa>=2){
	cout<<"from PndSttMvdTracking, evt. "<<IVOLTE<<", Mvd only track, Pivoting hit :\n";
	if(ListType[0]==0){ cout<<"\tPixel n. "<<List[0];}
	else if (ListType[0]==1){ cout<<"\tStrip n. "<<List[0];}
	cout<<", X "<<AloneX[0]<<", Y "<<AloneY[0]<<endl;
	cout<<"\tother hits :"<<endl;
	for(int h =1;h<nalone;h++){
	  if(ListType[h]==0){
		cout<<"\tPixel hit n. "<<List[h]<<", X "<<AloneX[h]<<", Y "<<AloneY[h]<<endl;
	  } else {
		cout<<"\tStrip hit n. "<<List[h]<<", X "<<AloneX[h]<<", Y "<<AloneY[h]<<endl;
	  }
	}
		cout<<" rotation angle prima = "<<rotationangle<<endl;
}
//--------------------


			// here at this point there are no Stt hits associated.
			RefitMvdStt(
				nalone,
				&fListTrackCandHit[nTotalCandidates][0],
				&fListTrackCandHitType[nTotalCandidates][0],
				info,
				rotationangle,
				trajectory_vertex,
				iexcl,

				&fALFA[nTotalCandidates], // output of the fit
				&fBETA[nTotalCandidates], // output of the fit
				&fGAMMA[nTotalCandidates],// output of the fit
				&status[nTotalCandidates]// fit status; true = success.
					);

				fOx[nTotalCandidates] = -fALFA[nTotalCandidates]/2.;
				fOy[nTotalCandidates] = -fBETA[nTotalCandidates]/2.;
				fR[nTotalCandidates] =
					fOx[nTotalCandidates]*fOx[nTotalCandidates]
					+fOy[nTotalCandidates]*fOy[nTotalCandidates]
					-fGAMMA[nTotalCandidates];
				if( fR[nTotalCandidates] > 0. ) fR[nTotalCandidates]=
					sqrt(fR[nTotalCandidates]) ;
				else  fR[nTotalCandidates]=0.;
			if(status[nTotalCandidates]&& fR[nTotalCandidates] > RSTRAWDETECTORMIN/2.){
				FI0[nTotalCandidates] =
				    atan2(-fOy[nTotalCandidates], -fOx[nTotalCandidates]);
				if( FI0[nTotalCandidates] < 0. )
					FI0[nTotalCandidates]+= 2.*PI;

				fnSttParHitsinTrack[nTotalCandidates]=0;
				fnSttSkewHitsinTrack[nTotalCandidates]=0;
				fnMvdPixelHitsinTrack[nTotalCandidates]=0;
				fnMvdStripHitsinTrack[nTotalCandidates]=0;
				for(j=0; j<nalone;j++){
					if(ListType[j]==0){
						fListMvdPixelHitsinTrack[nTotalCandidates][
							fnMvdPixelHitsinTrack[nTotalCandidates]
										]=List[j];
						fnMvdPixelHitsinTrack[nTotalCandidates]++;
					}else{
						fListMvdStripHitsinTrack[nTotalCandidates][
							fnMvdStripHitsinTrack[nTotalCandidates]
										]=List[j];
						fnMvdStripHitsinTrack[nTotalCandidates]++;
					}
				}  //   end of  for(j=1; j<nalone;j++)

				fnTrackCandHit[nTotalCandidates]=nalone;

				keepit[nTotalCandidates]=true;

				Mvdhits[nTotalCandidates]=true;
				nTotalCandidates++;
				nRemainingCandidates++;

			} // end of  if(status[nTotalCandidates])

		} // end of  if( nalone == 3 )

	} //end of  for(i=0; i<fnMvdTrackCand; i++)


//------------------ trying to attach a SciTil hit to the track.

	for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++){
//		if(!keepit[ncand]) continue;
		fnSciTilHitsinTrack[ncand]=0;	/// for the time being.
	}

//------------------ end trying to attach a SciTil hit to the track.

//     ordering all the hits belonging to the new candidate tracks, by increasing fR;
//     loading    fListTrackCandHit. The ordering is necessary here because the charge
//     has to be calculated.



	Ordering_Loading_ListTrackCandHit(
		keepit,
		nSttTrackCand,
		nTotalCandidates,
		info,
		Trajectory_Start,
		Charge,
		SchosenSkew  // here this array is irrelevant since there are not yet
			// any Stt hits in Track Cand.
		);

//-------------------- end of ordering

// -----------------  finding the charge of those track found starting from Mvd

 int dim = MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK+
			MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK;

 int nMvdOnly;
 Double_t
	esse[dim],
	Sini,
	Slast;


	for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++){
	    if(!keepit[ncand]) continue;

	// S[]  is used later in the SZ fit.

	    for(i=0, nMvdOnly=0; i<fnTrackCandHit[ncand]; i++){

		if( fListTrackCandHitType[ncand][i] == 0){ //  Pixel
			S[nMvdOnly] = atan2(
					fYMvdPixel[fListTrackCandHit[ncand][i]]-fOy[ncand],
					fXMvdPixel[fListTrackCandHit[ncand][i]]-fOx[ncand]
					);
			if(S[nMvdOnly]<0.) S[nMvdOnly] +=2.*PI;
			if(S[nMvdOnly]<0.) S[nMvdOnly] =0.;
			esse[i]=S[nMvdOnly];
			ZED[nMvdOnly] = fZMvdPixel[fListTrackCandHit[ncand][i]];
			// DriftRadius is set conventionally at -1, for later use in the SZ fit;
			// the error on the point used in the fit is ErrorDriftRadius and this
			// is overestimated to be  1cm.
			DriftRadius[nMvdOnly]= -1. ;
			ErrorDriftRadius[nMvdOnly]= 1. ;
			nMvdOnly++;
		} else if (fListTrackCandHitType[ncand][i] == 1) {  //  Strip
			S[nMvdOnly] = atan2(
					fYMvdStrip[fListTrackCandHit[ncand][i]]-fOy[ncand],
					fXMvdStrip[fListTrackCandHit[ncand][i]]-fOx[ncand]
					);
			if(S[nMvdOnly]<0.) S[nMvdOnly] +=2.*PI;
			if(S[nMvdOnly]<0.) S[nMvdOnly] =0.;
			esse[i]=S[nMvdOnly];
			ZED[nMvdOnly] = fZMvdStrip[fListTrackCandHit[ncand][i]];
			// DriftRadius is set conventionally at -1, for later use in the SZ fit;
			// the error on the point used in the fit is ErrorDriftRadius and this
			// is overestimated to be  1cm.
			DriftRadius[nMvdOnly]= -1. ;
			ErrorDriftRadius[nMvdOnly]= 1. ;
			nMvdOnly++;
		}   //   end of    if( fListTrackCandHitType[ncand][0] == 0)


	    } // end of  for(i=0; i<fnTrackCandHit[ncand]; i++)

	//  fixing possible discontinuity between fi<2*PI and fi>0.

	    flag = false;
	    for(i=1;i<fnTrackCandHit[ncand]; i++){
		if(fabs(esse[i-1] - esse[i]) > 1.5*PI ){
			flag=true;
			break;
		}
	    }
	    if( flag ){
		if( esse[0] < PI) Sini = esse[0]+2.*PI;
		if( esse[fnTrackCandHit[ncand]-1] < PI) Slast=esse[fnTrackCandHit[ncand]-1]+2.*PI;
	    } else {
		Sini = esse[0];
		Slast= esse[fnTrackCandHit[ncand]-1];
	    }

	    if( Sini > Slast) {
		Charge[ncand] =  1;
	    } else {
		Charge[ncand] =  -1;
	    }



	// it is necessary first to calculate Fi_low_limit and Fi_up_limit.

	FindingParallelTrackAngularRange(
		fOx[ncand],
		fOy[ncand],
		fR[ncand],
		Charge[ncand],
		&Fi_low_limit[ncand],	// Fi (in XY Helix frame) lower limit using
					// the Stt detector minimum/maximum radius
					// Fi_low_limit is ALWAYS between 0. and 2PI
		&Fi_up_limit[ncand],	// Fi (in XY Helix frame) upper limit using
					// the Stt detector maximum/minimum radius
					// Fi_up_limit is ALWAYS > Fi_low_limit and
					// possibly > 2PI.
		&statusflag[ncand],// vector; =0, all well; =1, track contained completely between RMin
			   // and RMax; = -1 track contained within RMin; =-2 track outside RMax.
		RSTRAWDETECTORMIN,
		RSTRAWDETECTORMAX
						);
		// discard tracks completely outside RMin (since they are formed by Mvd hits).
		if(statusflag[ncand]== -2 || statusflag[ncand]== 1){
			keepit[ncand]=false;
			continue;
		}
//  take care of possible discontinuities at 0.

	FixDiscontinuitiesFiangleinSZplane(
		fnTrackCandHit[ncand],	//  n hits in track. 
		S,	// S can be modified by +-2*PI if necessary.
		&FI0[ncand],	// this remains unchanged.
		Charge[ncand]	// this remains unchanged.
				);




// ----------------------------- fit in SZ with the Mvd tracks
		resultFitSZagain[ncand] = fit.FitSZspace(
				nMvdOnly,	// n. hits to be fitted
				tS,
				tZED,
				tDriftRadius,
				tErrorDriftRadius,
				FI0[ncand],
				MAXSKEWHITSINFIT,	// maximum STT Skew hits in fit;
				&emme,
				0   // IVOLTE
						);
		if( resultFitSZagain[ncand]==1){
			KAPPA[ncand] = emme;
			GoodSkewFit[ncand] = true;
		} else {
			keepit[ncand]=false;
			GoodSkewFit[ncand] = false;
		}

	} //   end of  for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++)

// ----------------------------- end fit in SZ with the Mvd tracks




//---  try to associate of parallel Stt  straw  hits to this new tracks.


	CollectParSttHitsagain(
			keepit,
			Mvdhits,
			info,
			nSttParHit,
			nSttTrackCand,	// starting candidate n. (included).
			nTotalCandidates, // ending candidate n. (excluded).
			KAPPA,
			FI0,
			Fi_low_limit,
			Fi_up_limit,
			fnSttParHitsinTrack, // input and output
			fListSttParHitsinTrack // input and output
			);
	// limit the total # Stt hits to MAXSTTHITSINTRACK
   for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++){
	if(!keepit[ncand]) continue;
	if( fnSttParHitsinTrack[ncand] > MAXSTTHITSINTRACK ) {
	  fnSttParHitsinTrack[ncand]=MAXSTTHITSINTRACK;
	}
   } // end of   for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++)


   //     try to attach skew hits to the new tracks (the result can also be 0).
   for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++){
	if(!keepit[ncand]) continue;
	// when statusflag[ncand]<0 the track does not intersect Stt region.
	if( statusflag[ncand]<0) continue;

	// the arrays  Temporary....   are related to the current candidate; they loose meaning
	// as soon as the loop over the candidate tracks finishes;

	   fnSttSkewHitsinTrack[ncand] = AssociateSkewHitsToXYTrack(
		fInclusionListStt,
		nSttSkewHit,
		fListSttSkewHits,
		fOx[ncand],   //  input : X of center of XY plane circle
		fOy[ncand],   //  input : Y of center of XY plane circle
		fR[ncand],   //  input : Radius of XY plane circle
		info,
		WDX,
		WDY,
		WDZ,
		Fi_low_limit[ncand],	// Fi (in Helix XY frame) lower limit using the Stt detector minimum/maximum radius
		Fi_up_limit[ncand],	// Fi (in Helix XY frame) upper limit using the Stt detector maximum/minimum radius
		Charge[ncand],
		TemporarySkewList, // output : TemporarySkewList[*][0] =  list of selected skew hits
				   // numbers(in skew numbering);
				   // TemporarySkewList[*][1] =  solution number (0 or 1) as per
				   // the calculateintersections  method;
	// the info on TemporaryS, TemporaryZ, TemporaryZDrift, TemporaryZErrorafterTilt  as follows :
	// if   i  is the (original) Skew Hit number, and ii (0 or 1) is
	// the solution, then the infos are stored in the   i + ii*MAXSTTHITS location;
		TemporaryS,       //  output,  S coordinate of selected Skew hit
		TemporaryZ,       //  output,  Z coordinate of selected Skew hit (center wire)
		TemporaryZDrift,   //  output,  drift distance IN Z DIRECTION only, of selected Skew hit
		TemporaryZErrorafterTilt   //  output,  Radius taking into account the tilt, IN Z DIRECTION only, of selected Skew hit
								);
	   // limit the total # Stt hits to MAXSTTHITSINTRACK
	   if( fnSttSkewHitsinTrack[ncand]+fnSttParHitsinTrack[ncand] > MAXSTTHITSINTRACK ) {
		if(MAXSTTHITSINTRACK > fnSttParHitsinTrack[ncand])
		  fnSttSkewHitsinTrack[ncand]=MAXSTTHITSINTRACK-fnSttParHitsinTrack[ncand];
		else fnSttSkewHitsinTrack[ncand]=0;
	   }

	   for(j=0;j<fnSttSkewHitsinTrack[ncand];j++)
	   {
		fListSttSkewHitsinTrack[ncand][j]=TemporarySkewList[j][0];
		fListSttSkewHitsinTrackSolution[ncand][j]=TemporarySkewList[j][1];
		SchosenSkew[ncand][fListSttSkewHitsinTrack[ncand][j]] = TemporaryS[j];
	   }


  //  ------------------- end of attachment of skew hits to the new tracks.



   } // end of   for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++)



//     ordering again all the hits belonging to the new candidate tracks, by increasing fR;
//     loading    fListTrackCandHit.

	Ordering_Loading_ListTrackCandHit(
		keepit,
		nSttTrackCand,
		nTotalCandidates,
		info,
		Trajectory_Start,
		Charge,
		SchosenSkew
		);

//-------------------- end of ordering

	// adding at the end the SciTil hit (if present).

	for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++){
		if(!keepit[ncand]) continue;
		if( fnSciTilHitsinTrack[ncand]==1) {
			i=fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+
			  fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand];
			fListTrackCandHit[ncand][i] = fListSciTilHitsinTrack[ncand][0];
			fListTrackCandHitType[ncand][i] = 1001;
		}
	}  // end of for(ncand=0; ncand< nTotalCandidates; ncand++)


//-----------  cleanup of the new tracks found



  if(fYesClean){
     for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++){
	if(!keepit[ncand]) continue;

	if ( !Cleaner.TrackCleanup(
		APOTEMAMAXINNERPARSTRAW,
		APOTEMAMAXSKEWSTRAW,
		APOTEMAMINOUTERPARSTRAW,
		APOTEMAMINSKEWSTRAW,
		&SchosenSkew[ncand][0],
		Charge[ncand],
		FI0[ncand],
		gap,
		info,
		istampa,
		IVOLTE,
		KAPPA[ncand],
		&fListSttParHitsinTrack[ncand][0],
		&fListSttSkewHitsinTrack[ncand][0],
		MAXSTTHITS,
		fnSttParHitsinTrack[ncand],
		fnSttSkewHitsinTrack[ncand],
		fOx[ncand],
		fOy[ncand],
		fR[ncand],
		RSTRAWDETECTORMAX,
		RSTRAWDETECTORMIN,
		SEMILENGTH_STRAIGHT,
		Start,
		STRAWRADIUS,
		ZCENTER_STRAIGHT
				) ) {
		keepit[ncand]=false;
	}  // end of if ( !TrackCleanup


     } // end of   for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++)
  }  // end of if(fYesClean)

//-----------  end of cleanup of new tracks found.


	} // end of  if(fMvdAloneTracking && fnMvdPixelHit+fnMvdStripHit>0)



// ----------------------------- end of finding Tracks starting with the Mvd tracks
// -------------------------------------------------------------------------------------
// -------------------------------------------------------------------------------------
// -------------------------------------------------------------------------------------

//-----------------------
if(istampa>1){
 for(ncand=0; ncand< nTotalCandidates; ncand++) {
	cout<<"from PndTrkTracking.cxx, just before  LoadPndTrack_TrackCand, cand. # "
	<<ncand<<", keepit = "<<keepit[ncand];
	if(keepit[ncand]){
		  cout<<", # Pixels in cand = "<<fnMvdPixelHitsinTrack[ncand]<<" e loro lista :\n";
		for(int gg=0; gg<fnMvdPixelHitsinTrack[ncand]; gg++){
		  cout<<"\tPixel Hit # "<<fListMvdPixelHitsinTrack[ncand][gg]<<endl;
		}
		cout<<ncand<<", # Strips in cand = "<<fnMvdStripHitsinTrack[ncand]<<" e loro lista :\n";
        	for(int gg=0; gg<fnMvdStripHitsinTrack[ncand]; gg++){
                	cout<<"\tPixel Hit # "<<fListMvdStripHitsinTrack[ncand][gg]<<endl;
        	}
	} else {
		cout<<";\n";
	}  // end of  if(keepit[ncand])

 }
}


        for(ncand=0; ncand< nTotalCandidates; ncand++){
                if(!keepit[ncand]) continue;
                for(int jj=0; jj<fnTrackCandHit[ncand]; jj++){
                        switch (fListTrackCandHitType[ncand][jj]){
                                case 0:
                                        if ( fListTrackCandHit[ncand][jj] >=fnMvdPixelHit ){
                                                cout<<"from PndTrkTracking, problem event "<<IVOLTE<<", cand. "<<ncand<<", hit n. "<<
                                                fListTrackCandHit[ncand][jj]<<"  type Mvd Pixel; (n. Mvd Pixel Hit = "<<fnMvdPixelHit<<" );\n";
                                        }
                                        break;

                                case 1:
                                        if ( fListTrackCandHit[ncand][jj] >=fnMvdStripHit ){
                                                cout<<"from PndTrkTracking, problem event "<<IVOLTE<<", cand. "<<ncand<<", hit n. "<<
                                                fListTrackCandHit[ncand][jj]<<"  type Mvd Strip; (n. Mvd Strip Hit = "<<fnMvdStripHit  <<" );\n";
                                        }
                                        break;

                                case 2:
                                        if ( fListTrackCandHit[ncand][jj] >=nSttHit ){
                                                cout<<"from PndTrkTracking, problem event "<<IVOLTE<<", cand. "<<ncand<<", hit n. "<<
                                                fListTrackCandHit[ncand][jj]<<"  type Stt || hit; (nSttHit = "<<nSttHit<<" );\n";
                                        }
                                        break;
                                case 3:
                                        if ( fListTrackCandHit[ncand][jj] >=nSttHit ){
                                                cout<<"from PndTrkTracking, problem event "<<IVOLTE<<", cand. "<<ncand<<", hit n. "<<
                                                fListTrackCandHit[ncand][jj]<<"  type Stt // hit; (nSttHit = "<<nSttHit<<" );\n";
                                        }
                                        break;
                                case 1001:
                                        if ( fListTrackCandHit[ncand][jj] >=fnSciTilHits ){
                                                cout<<"from PndTrkTracking, problem event "<<IVOLTE<<", cand. "<<ncand<<", hit n. "<<
                                                fListTrackCandHit[ncand][jj]<<"  type SciTil hit; (n. SciTi Hit = "<<fnSciTilHits<<" );\n";
                                        }
                                        break;

                        }
                }
        } // end of   for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++)

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















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

//-------  load the new PndTrackCand ; each track has the STT and the Mvd hits associated
//-------  also load the new PndTrack ; each track has the STT and the Mvd hits associated

 LoadPndTrack_TrackCand(
	keepit,
	SttSZfit,
	nTotalCandidates,
	Charge,
	nSttTrackCand,
	FI0,
	KAPPA,
	info,
	SchosenSkew,
	ZchosenSkew
	);


//------------ section with comparison MC Mvd hits - associated hits to a certain track




 int dim1,  dim2  , dim3 , dim4, dim5, dim6, dim7;

 if(nSttHit == 0 ) { dim1=1; } else { dim1 = MAXTRACKSPEREVENT*nSttHit; };

 Short_t
	MCParalAloneList[dim1],
	MCSkewAloneList[dim1];

 if(nTotalCandidates == 0 ) {
	dim1 = 1;
	dim2 = 1;
	dim3 = 1;
	dim4 = 1;
	dim5 = 1;
	dim6 = 1;
	dim7 = 1;
 } else {
	dim1 = nTotalCandidates;
	dim2 = nTotalCandidates*MAXMVDPIXELHITSINTRACK;
	dim3 = nTotalCandidates*MAXMVDSTRIPHITSINTRACK;
	dim4 = nTotalCandidates*fnMvdPixelHit;
	dim5 = nTotalCandidates*fnMvdStripHit;
	dim6 = nTotalCandidates*MAXSCITILHITSINTRACK;
	dim7 = nTotalCandidates*fnSciTilHits;
 }

 Short_t
	nMvdPixelCommon[dim1],
	MvdPixelCommonList[dim2],
	nMvdPixelSpuriinTrack[dim1],
	MvdPixelSpuriList[dim2],
	nMCMvdPixelAlone[dim1],
	MCMvdPixelAloneList[dim4],

	nMvdStripCommon[dim1],
	MvdStripCommonList[dim3],
	nMvdStripSpuriinTrack[dim1],
	MvdStripSpuriList[dim3],
	nMCMvdStripAlone[dim1],
	MCMvdStripAloneList[dim5],

	nSciTilCommon[dim1],
	SciTilCommonList[dim6],
	nSciTilSpuriinTrack[dim1],
	SciTilSpuriList[dim6],
	nMCSciTilAlone[dim1],
	MCSciTilAloneList[dim7];


 if( doMcComparison && nTotalCandidates > 0){
	// make the struct for the data to pass to the
	// method PndTrkComparisonMCtruth::ComparisonwithMC ;
	PndTrkComparisonMCtruth_io_Data ioData;
	// load the structure;

	ioData.Bfield = BFIELD;
	ioData.Charge = tCharge;
	ioData.Cvel = CVEL;
	ioData.daTrackFoundaTrackMC = tdaTrackFoundaTrackMC;
	ioData.DIMENSIONSciTil = DIMENSIONSCITIL;
	ioData.Errorsqpixel = ERRORSQPIXEL;
	ioData.Errorsqstrip = ERRORSQSTRIP;
	ioData.FI0 = FI0;
	ioData.fMCTrackArray = fMCTrackArray;
	ioData.fMvdMCPointArray = fMvdMCPointArray;
	ioData.fSciTilMaxNumber = fSciTilMaxNumber;
	ioData.fSciTHitArray = fSciTHitArray;
	ioData.fSciTPointArray = fSciTPointArray;
	ioData.fSttPointArray = fSttPointArray;
	ioData.HANDLE = HANDLE,
	ioData.HANDLE2 = HANDLE2,
	ioData.info = &info[0][0];
	ioData.istampa = istampa;
	ioData.IVOLTE = IVOLTE;
	ioData.KAPPA = KAPPA;
	ioData.keepit = tkeepit;
	ioData.InclusionListStt = fInclusionListStt;
	ioData.ListMvdPixelHitsinTrack = &fListMvdPixelHitsinTrack[0][0];
	ioData.ListMvdStripHitsinTrack = &fListMvdStripHitsinTrack[0][0];
	ioData.ListSciTilHitsinTrack = &fListSciTilHitsinTrack[0][0];
	ioData.ListSttParHitsinTrack = &fListSttParHitsinTrack[0][0];
	ioData.ListSttSkewHitsinTrack = &fListSttSkewHitsinTrack[0][0];
	ioData.ListTrackCandHit = &fListTrackCandHit[0][0];
	ioData.ListTrackCandHitType = &fListTrackCandHitType[0][0];
	ioData.MAXMCTRACKS = MAXMCTRACKS;
	ioData.MAXMVDPIXELHITS = MAXMVDPIXELHITS;
	ioData.MAXMVDPIXELHITSINTRACK = MAXMVDPIXELHITSINTRACK;
	ioData.Maxmvdmcpoints = MAXMVDMCPOINTS;
	ioData.MAXMVDSTRIPHITS = MAXMVDSTRIPHITS;
	ioData.MAXMVDSTRIPHITSINTRACK = MAXMVDSTRIPHITSINTRACK;
	ioData.MAXSCITILHITS = MAXSCITILHITS ;
	ioData.MAXSCITILHITSINTRACK = MAXSCITILHITSINTRACK;
	ioData.MAXSTTHITS = MAXSTTHITS;
	ioData.MAXSTTHITSINTRACK = MAXSTTHITSINTRACK;
	ioData.MAXTRACKSPEREVENT = MAXTRACKSPEREVENT;
	ioData.MCMvdPixelAloneList = MCMvdPixelAloneList;
	ioData.MCMvdStripAloneList = MCMvdStripAloneList;
	ioData.MCParalAloneList = MCParalAloneList;
	ioData.MCSciTilAloneList = MCSciTilAloneList;
	ioData.MCSkewAloneList = MCSkewAloneList;
	ioData.MCSkewAloneX = fMCSkewAloneX ;
	ioData.MCSkewAloneY = fMCSkewAloneY ;
	ioData.MvdPixelCommonList = MvdPixelCommonList;
	ioData.MvdPixelSpuriList = MvdPixelSpuriList;
	ioData.MvdStripCommonList = MvdStripCommonList;
	ioData.MvdStripSpuriList = MvdStripSpuriList;
	ioData.nHitsInMCTrack = tnHitsInMCTrack;
	ioData.nHitsInSciTile = nHitsInSciTile;
	ioData.nMCMvdPixelAlone = nMCMvdPixelAlone;
	ioData.nMCMvdStripAlone = nMCMvdStripAlone;
	ioData.nMCParalAlone = tnMCParalAlone;
	ioData.nMCSciTilAlone = nMCSciTilAlone;
	ioData.nMCSkewAlone = tnMCSkewAlone;
	ioData.nMvdPixelCommon = nMvdPixelCommon;
	ioData.nMvdPixelHitsinTrack = fnMvdPixelHitsinTrack;
	ioData.nMvdStripHitsinTrack = fnMvdStripHitsinTrack;
	ioData.nMvdPixelHit = fnMvdPixelHit;
	ioData.nMvdPixelSpuriinTrack = nMvdPixelSpuriinTrack;
	ioData.nMvdStripCommon = nMvdStripCommon;
	ioData.nMvdStripHit = fnMvdStripHit;
	ioData.nMvdStripSpuriinTrack = nMvdStripSpuriinTrack;
	ioData.nParalCommon = tnParalCommon;
	ioData.nSciTilCommon = nSciTilCommon;
	ioData.nSciTilHits = fnSciTilHits;
	ioData.nSciTilHitsinTrack = fnSciTilHitsinTrack ;
	ioData.nSciTilSpuriinTrack = nSciTilSpuriinTrack ;
	ioData.nSkewCommon = tnSkewCommon;
	ioData.nSkewHitsInMCTrack = tnSkewHitsInMCTrack;
	ioData.nSpuriParinTrack = tnSpuriParinTrack;
	ioData.nSpuriSkewinTrack = tnSpuriSkewinTrack;
	ioData.nSttHit = nSttHit;
	ioData.nSttParHitsinTrack = fnSttParHitsinTrack;
	ioData.nSttSkewHitsinTrack = fnSttSkewHitsinTrack;
	ioData.nTotalCandidates = nTotalCandidates;
	ioData.OriginalSciTilList = &OriginalSciTilList[0][0];
	ioData.Ox = fOx;
	ioData.Oy = fOy;
	ioData.ParalCommonList = tParalCommonList;
	ioData.ParSpuriList = tParSpuriList;
	ioData.R = fR;
	ioData.refindexMvdPixel = frefindexMvdPixel;
	ioData.refindexMvdStrip = frefindexMvdStrip;
	ioData.resultFitSZagain = tresultFitSZagain;
	ioData.SciTilCommonList = SciTilCommonList;
	ioData.SciTilSpuriList = SciTilSpuriList;
	ioData.SkewCommonList = tSkewCommonList;
	ioData.SkewSpuriList = tSkewSpuriList;
	ioData.SttSZfit = tSttSZfit;
	ioData.XMvdPixel = fXMvdPixel;
	ioData.XMvdStrip = fXMvdStrip;
	ioData.XSciTilCenter = fpSciTilx;
	ioData.YMvdPixel = fYMvdPixel;
	ioData.YMvdStrip = fYMvdStrip;
	ioData.YSciTilCenter = fpSciTily;
	ioData.ZMvdPixel = fZMvdPixel;
	ioData.ZMvdStrip = fZMvdStrip;
	ioData.ZSciTilCenter = fpSciTilz;

	// class for the MC comparison;
	PndTrkComparisonMCtruth cmp;
	fnMCTracks = cmp.ComparisonwithMC( ioData);

 }

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

 // write the Macro for visualization of tracks and hits;

 if(iplotta && IVOLTE< fNevents_to_plot){

	// the following initialization is necessary when the MC comparison
	// is not done just above (when doMcComparison=false). In this case
	// in fact it is necessary to have  the arrays nParalCommon, nSpuriParinTrack  etc.etc.
	// set at 0  otherwise some WriteMacro  methods  crash;
	if(!doMcComparison){
	 for(i=0; i<nTotalCandidates;i++){
	  nParalCommon[i]=0; nSpuriParinTrack[i]=0; nMCParalAlone[i]=0;
	  nSkewCommon[i]=0; nSpuriSkewinTrack[i]=0; nMCSkewAlone[i]=0;
	  nMvdPixelCommon[i]=0; nMvdPixelSpuriinTrack[i]=0; nMCMvdPixelAlone[i]=0;
	  nMvdStripCommon[i]=0; nMvdStripSpuriinTrack[i]=0; nMCMvdStripAlone[i]=0;
	  nSciTilCommon[i]=0; nSciTilSpuriinTrack[i]=0; nMCSciTilAlone[i]=0;

	 }  // end for(i=0; i<nTotalCandidates;i++)

	}


	PndTrkPlotMacros mymacro;
	PndTrkPlotMacros_InputData In_Put;
	In_Put.apotemamaxinnerparstraw = APOTEMAMAXINNERPARSTRAW ;
	In_Put.apotemamaxskewstraw = APOTEMAMAXSKEWSTRAW ;
	In_Put.apotemaminouterparstraw = APOTEMAMINOUTERPARSTRAW ;
	In_Put.apotemaminskewstraw = APOTEMAMINSKEWSTRAW ;
	In_Put.bfield = BFIELD ;
	In_Put.Charge = tCharge ;
	In_Put.cvel = CVEL ;
	In_Put.daTrackFoundaTrackMC = tdaTrackFoundaTrackMC ;
	In_Put.dimensionscitil = DIMENSIONSCITIL ;
	In_Put.doMcComparison = doMcComparison ;
	In_Put.FI0 = FI0 ;
	In_Put.fMCTrackArray = fMCTrackArray ;
	In_Put.fSttPointArray = fSttPointArray ;
	In_Put.info =  &info[0][0] ;
	In_Put.IVOLTE = IVOLTE ;
	In_Put.KAPPA = KAPPA ;
	In_Put.keepit = tkeepit ;
	In_Put.InclusionListSciTil = fInclusionListSciTil ;
	In_Put.InclusionListStt = SingleHitListStt ;

	In_Put.istampa = istampa ;
	In_Put.ListMvdPixelHitsinTrack = &fListMvdPixelHitsinTrack[0][0] ;
	In_Put.ListMvdStripHitsinTrack = &fListMvdStripHitsinTrack[0][0] ;
	In_Put.ListSciTilHitsinTrack = &fListSciTilHitsinTrack[0][0] ;
	In_Put.ListSttParHitsinTrack = &fListSttParHitsinTrack[0][0] ;
	In_Put.ListSttSkewHitsinTrack = &fListSttSkewHitsinTrack[0][0] ;
	In_Put.ListTrackCandHit = &fListTrackCandHit[0][0] ;
	In_Put.ListTrackCandHitType = &fListTrackCandHitType[0][0] ;
	In_Put.MAXMCTRACKS = MAXMCTRACKS ;
	In_Put.MAXMVDPIXELHITS = MAXMVDPIXELHITS ;
	In_Put.MAXMVDPIXELHITSINTRACK = MAXMVDPIXELHITSINTRACK ;
	In_Put.MAXMVDSTRIPHITS = MAXMVDSTRIPHITS ;
	In_Put.MAXMVDSTRIPHITSINTRACK = MAXMVDSTRIPHITSINTRACK ;
	In_Put.MAXSCITILHITSINTRACK = MAXSCITILHITSINTRACK ;
	In_Put.MAXSCITILHITS = MAXSCITILHITS ;
	In_Put.MAXSTTHITS = MAXSTTHITS ;
	In_Put.MAXSTTHITSINTRACK = MAXSTTHITSINTRACK ;
	In_Put.MAXTRACKSPEREVENT = MAXTRACKSPEREVENT ;
	In_Put.MCMvdPixelAloneList = MCMvdPixelAloneList ;
	In_Put.MCMvdStripAloneList = MCMvdStripAloneList ;
	In_Put.MCParalAloneList = MCParalAloneList ;
	In_Put.MCSciTilAloneList = MCSciTilAloneList;
	In_Put.MCSkewAloneList = MCSkewAloneList ;
	In_Put.MCSkewAloneX = fMCSkewAloneX ;
	In_Put.MCSkewAloneY = fMCSkewAloneY ;
	In_Put.MvdPixelCommonList = MvdPixelCommonList ;
	In_Put.MvdPixelSpuriList = MvdPixelSpuriList ;
	In_Put.MvdStripCommonList = MvdStripCommonList ;
	In_Put.MvdStripSpuriList = MvdStripSpuriList ;
	In_Put.NFIDIVCONFORMAL =  NFIDIVCONFORMAL;
	In_Put.nMCMvdPixelAlone = nMCMvdPixelAlone ;
	In_Put.nMCMvdStripAlone = nMCMvdStripAlone ;
	In_Put.nMCParalAlone = tnMCParalAlone ;
	In_Put.nMCSciTilAlone = nMCSciTilAlone;
	In_Put.nMCSkewAlone = tnMCSkewAlone ;
	In_Put.nMCTracks = fnMCTracks ;
	In_Put.nMvdPixelCommon = nMvdPixelCommon ;
	In_Put.nMvdPixelHit = fnMvdPixelHit ;
	In_Put.nMvdPixelHitsinTrack = fnMvdPixelHitsinTrack ;
	In_Put.nMvdPixelSpuriinTrack = nMvdPixelSpuriinTrack ;
	In_Put.nMvdStripCommon = nMvdStripCommon ;
	In_Put.nMvdStripHit = fnMvdStripHit ;
	In_Put.nMvdStripHitsinTrack = fnMvdStripHitsinTrack ;
	In_Put.nMvdStripSpuriinTrack = nMvdStripSpuriinTrack ;
	In_Put.NRDIVCONFORMAL = NRDIVCONFORMAL;
	In_Put.nParalCommon = tnParalCommon ;
	In_Put.nSciTilCommon = nSciTilCommon ;
	In_Put.nSciTilHits = fnSciTilHits ;
	In_Put.nSciTilHitsinTrack = fnSciTilHitsinTrack ;
	In_Put.nSciTilSpuriinTrack =  nSciTilSpuriinTrack;
	In_Put.nSkewCommon = tnSkewCommon ;
	In_Put.nSpuriParinTrack = tnSpuriParinTrack ;
	In_Put.nSttHit = nSttHit ;
	In_Put.nSttParHit = nSttParHit ;
	In_Put.nSttParHitsinTrack = fnSttParHitsinTrack ;
	In_Put.nSttSkewHit = nSttSkewHit ;
	In_Put.nSttSkewHitsinTrack = fnSttSkewHitsinTrack ;
	In_Put.nTotalCandidates = nTotalCandidates ;
	In_Put.nTrackCandHit = fnTrackCandHit ;
	In_Put.Ox = fOx ;
	In_Put.Oy = fOy ;
	In_Put.ParalCommonList = tParalCommonList ;
	In_Put.ParSpuriList = tParSpuriList ;
	In_Put.posizSciTil = &fposizSciTil[0][0] ;
	In_Put.R = fR ;
	In_Put.radiaConf = fradiaConf;
	In_Put.rstrawdetectormax = RSTRAWDETECTORMAX ;
	In_Put.rstrawdetectormin = RSTRAWDETECTORMIN ;
	In_Put.SchosenSkew = &SchosenSkew[0][0] ;
	In_Put.SciTilCommonList = SciTilCommonList ;
	In_Put.SciTilSpuriList = SciTilSpuriList ;
	In_Put.sigmaXMvdPixel = fsigmaXMvdPixel ;
	In_Put.sigmaXMvdStrip = fsigmaXMvdStrip ;
	In_Put.sigmaYMvdPixel = fsigmaYMvdPixel ;
	In_Put.sigmaYMvdStrip = fsigmaYMvdStrip ;
	In_Put.SkewCommonList = tSkewCommonList ;
	In_Put.verticalgap = VERTICALGAP ;
	In_Put.XMvdPixel = fXMvdPixel ;
	In_Put.XMvdStrip = fXMvdStrip ;
	In_Put.YMvdPixel = fYMvdPixel ;
	In_Put.YMvdStrip = fYMvdStrip ;
	In_Put.WDX = WDX ;
	In_Put.WDY = WDY ;
	In_Put.WDZ = WDZ ;
	In_Put.ZMvdPixel = fZMvdPixel ;
	In_Put.ZMvdStrip = fZMvdStrip ;

	mymacro.WriteAllMacros(
	In_Put
				);
 }

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




 return;

}


//---------------------- end of   PndTrkTracking::Exec


//----------begin of function PndTrkTracking::AssociateBetterAfterFitSkewHitsToXYTrack

  Short_t PndTrkTracking::AssociateBetterAfterFitSkewHitsToXYTrack(
                   Short_t TemporarynSttSkewhitinTrack,  //  input
                   Short_t SkewList[][2], // input,  list of selected skew hits (in skew numbering)
                   Double_t *S,       //  input,  S coordinate of selected Skew hit
                   Double_t *Z,       //  input,  Z coordinate of selected Skew hit
                   Double_t *ZDrift,  //  input,  drift distance IN Z DIRECTION only, of selected Skew hit
                   Double_t *ZErrorafterTilt,   //  input,  Radius taking into account the tilt, IN Z DIRECTION only, of selected Skew hit
                   Double_t KAPPA,    // input, KAPPA result of fit
                   Double_t FI0,    // input, FI0 result of fit
                   Short_t *tempore,  //  output result, associated skew hits
                   Double_t *temporeS,  //  output, associated skew hit  S
                   Double_t *temporeZ,  //  output, associated skew hits Z
                   Double_t *temporeZDrift,  //  output, associated skew hit Z drift
                   Double_t *temporeZErrorafterTilt,  //  output, associated skew hits Z error after tilt
                   Short_t  *STATUS   // output
                                                     )
 {

    Short_t NAssociated;
    Short_t  sign;
    Int_t i, j, i1, ii, iii,  Kincl, nlow, nup;

    Double_t bbb,
             tempZ[2],
             zmin, zmax, deltaz,
		zdist[2],
		zdist1,
		zdist2;

    Double_t allowed_distance = 4.*STRAWRADIUS/sin(STRAW_SKEW_INCLINATION*PI/180.);

    if(fabs(KAPPA)<1.e-20) {
      *STATUS=-1;
      return 0;
    }

    NAssociated=0;
    if(KAPPA>0) {
     zmin = -FI0/KAPPA;
     zmax = (2.*PI-FI0)/KAPPA;
    }  else {
     zmax = -FI0/KAPPA;
     zmin = (2.*PI-FI0)/KAPPA;
    }
    deltaz = zmax-zmin;

     for(i=0; i<TemporarynSttSkewhitinTrack; i++){
       bbb=(S[i]-FI0)/KAPPA;
      for(sign=0;sign<=1; sign ++){
       tempZ[sign]=Z[i]+(2*sign-1)*ZDrift[i];
       if( tempZ[sign] > zmax ){
         tempZ[sign]=fmod( tempZ[sign]-zmax, deltaz) + zmin;
       } else if (tempZ[sign]<zmin){
         tempZ[sign]=fmod( tempZ[sign]-zmin, deltaz) + zmax;
       }
	zdist1 = fabs( bbb - tempZ[sign]);
	zdist2 = deltaz- zdist1;
	if(zdist2<0.) zdist2 = 0.;  // protect against rounding errors.
	zdist[sign] = zdist1 < zdist2 ? zdist1 : zdist2;
      }  //  end of for(sign=0;sign<=1; sign++)
	zdist1 = zdist[0] < zdist[1] ? zdist[0] : zdist[1];
	if(  zdist1 < allowed_distance ){
         tempore[NAssociated]=SkewList[ i ][0];
         temporeS[NAssociated]=S[i];
         temporeZ[NAssociated]=Z[i];
         temporeZDrift[NAssociated]=ZDrift[i];
         temporeZErrorafterTilt[NAssociated]=ZErrorafterTilt[i];
         NAssociated++;
       }
     }   //  end of for(i=0; i<TemporarynSttSkewhitinTrack; i++)

     *STATUS=0;
     return NAssociated;
 }

//----------end of function PndTrkTracking::AssociateBetterAfterFitSkewHitsToXYTrack



//----------begin of function PndTrkTracking::AssociateSkewHitsToXYTrack

Short_t PndTrkTracking::AssociateSkewHitsToXYTrack(
	bool *InclusionListSkew,
	Short_t NSkewhits,
	Short_t *infoskew,  // from skew numbering to original Stt hit numbering;
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t info[][7],
	Double_t *WDX,
	Double_t *WDY,
	Double_t *WDZ,
	Double_t Fi_low_limit,
	Double_t Fi_up_limit,
	Short_t  Charge,
	Short_t SkewList[][2], // output,list of selected skew hits (original numbering)
	Double_t *S,       //  output,  S coordinate of selected Skew hit
	Double_t *Z,       //  output,  Z coordinate of selected Skew hit
	Double_t *ZDrift,   //  output,  drift distance IN Z DIRECTION only,
			   // of selected Skew hit
	Double_t *ZErrorafterTilt   //  output,  Radius taking into account
				// the tilt, IN Z DIRECTION only, of selected Skew hit.
	)
 {

 int
	i,
	i1,
	ii,
	iii,
	j,
	Kincl,
	location,
	NAssociated,
	nlow,
	Nmax,
	Nmin,
	nup,
	STATUS;

 Double_t xmin , xmax, ymin, ymax,
           dx, dy, diff, dista, d1, d2,
           delta, deltax, deltay, deltaz, deltaS,
           factor,
           zmin, zmax, Smin, Smax, S1, S2,
           z1, z2, y1, y2,
           vx1, vy1, vz1, C0x1, C0y1, C0z1,
           aaa, bbb, ccc, angle, minor, major,
           distance, Rx, Ry, LL,
           Aellipsis1, Bellipsis1,fi1,
           fmin, fmax, offset, step,
           SkewInclWithRespectToS,
	Tilted_Radius,
	zpos, zpos1, zpos2,
           Tiltdirection1[2],
           zl[200],zu[200],
           POINTS1[6];


 PndTrkCTGeometryCalculations GeomCalculator;


      NAssociated=0;

       for( iii=0; iii< NSkewhits; iii++) {
         i = infoskew[iii];  // i  is the Hit number in original Skew Hit numbering;
         if( !InclusionListSkew[i]) continue;
         aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+WDZ[i]*WDZ[i]);
         vx1 = WDX[i]/aaa;
         vy1 = WDY[i]/aaa;
         vz1 = WDZ[i]/aaa;
         C0x1 = info[i][0];
         C0y1 = info[i][1];
         C0z1 = info[i][2];
       GeomCalculator.calculateintersections(Oxx,Oyy,Rr,C0x1,C0y1,C0z1,info[i][3],
                              vx1,vy1,vz1,
                              &STATUS,POINTS1);

       if(STATUS < 0 ) continue ;


       for( ii=0; ii<2; ii++){

	// ii =0, 1 is the solution number according to the algoritm of
	// the method   calculateintersections;
        j=3*ii;
        distance = sqrt(
                  (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + 
                  (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + 
                  (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) 
                            );


        Rx = POINTS1[j]-Oxx ;   //  x component Radial vector of cylinder of trajectory
        Ry = POINTS1[1+j]-Oyy ;   //  y direction Radial vector of cylinder of trajectory

        aaa = sqrt(Rx*Rx+Ry*Ry);
        SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ;
        SkewInclWithRespectToS /= Rr;
        bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1);
        //  the tilt direction of this ellipse is (1,0)  when major axis along Z direction
        if( bbb > 1.e-10){
           Tiltdirection1[0] = vz1/bbb;
           Tiltdirection1[1] = SkewInclWithRespectToS/bbb;
        } else {
           Tiltdirection1[0] = 1.;
           Tiltdirection1[1] = 0.;
        }

        LL = fabs(vx1*Rx + vy1*Ry);
        if( LL < 1.e-10) continue;

        Aellipsis1 = info[i][3]*aaa/LL;

        Bellipsis1 = info[i][3]/Rr;

	Tilted_Radius = STRAWRADIUS*aaa/LL;
	//   info[i][4]  =  semilength of the wire;

	// the following statement adds some margin to take into account the imprecision in the determination
	// of the Helix radius;
	dista = 2.*Aellipsis1 > Tilted_Radius ?  Tilted_Radius :  2.*Aellipsis1;

        if( distance >= info[i][4]+ dista ){
		 continue;
	}

	// store S, Z, ZDrift, ZErrorafterTilt  as follows :
	// if   i  is the (original) Skew Hit number, and ii (0 or 1) is
	// the solution, then the infos are stored in the   i + ii*MAXSTTHITS location;

	location = i + ii*MAXSTTHITS;


        S[location] = atan2(POINTS1[j+1]-Oyy, POINTS1[j]-Oxx) ;  // atan2 returns radians in (-pi and +pi]
        if( S[location] < 0.) S[location] += 2.*PI;


        if(  S[location] < Fi_low_limit) {
           if(  S[location]+2.*PI > Fi_up_limit)  continue;
        }  else if(  S[location] > Fi_up_limit) {
	   if(  S[location]- 2.*PI < Fi_low_limit)  continue;
        }



//---------------------------   end check


	Z[location] = POINTS1[j+2];
	ZDrift[location] = Aellipsis1*Tiltdirection1[0];
	ZErrorafterTilt[location] = 0.02*aaa*Tiltdirection1[0]/LL;

	// the list of the selected Skew skew hits instead, is stored sequentially;

	SkewList[NAssociated][0] = i;  // n. skew hit in original hit numbering
	SkewList[NAssociated][1] = ii;  //  solution 0 or solution 1 were accepted


// check if this skew hit doesn't "push out"  the most external parallel hit
// (see Gianluigi's logbook on page 251)

/*
        Double_t Zh1 = Z[location] - ZDrift[location];
        Double_t Zh2 = Z[location] + ZDrift[location];
        Double_t Sh1 = S[location] - Aellipsis1*Tiltdirection1[1];
        Double_t Sh2 = S[location] + Aellipsis1*Tiltdirection1[1];
        Double_t Zlast1 = (Fi_final_helix_referenceframe-Fi_initial_helix_referenceframe)*Zh1
                                  /(Sh1-Fi_initial_helix_referenceframe);
        Double_t Zlast2 = (Fi_final_helix_referenceframe-Fi_initial_helix_referenceframe)*Zh2
                                  /(Sh2-Fi_initial_helix_referenceframe);



        if( fabs(Zlast1 - ZCENTER_STRAIGHT) > SEMILENGTH_STRAIGHT
                                  &&
            fabs(Zlast2 - ZCENTER_STRAIGHT) > SEMILENGTH_STRAIGHT
           )   continue;
*/
        NAssociated++;

   }    //  end of    for( ii=0; ii<2; ii++)

  }   //   for( iii=0; iii< NSkewhits; iii++)

  return NAssociated;



}

//----------end of function PndTrkTracking::AssociateSkewHitsToXYTrack


//------------------------- begin of function  PndTrkTracking::CollectParSttHitsagain

 void PndTrkTracking::CollectParSttHitsagain(
	Vec <bool>& keepit,
	Vec <bool>& Mvdhits,
	Double_t info[][7],
	Short_t nSttParHit,
	// starting investigation from candidate n. StartTrackCand
	Short_t StartTrackCand,
	//  until candidate n. EndTrackCand-1 included.
	Short_t EndTrackCand,
	Double_t *KAPPA,
	Double_t *FI0,
	Double_t *Fi_low_limit,
	Double_t *Fi_up_limit,
	Short_t *nParHitsinTrack, // input/output
	Short_t ListParHitsinTrack[][MAXSTTHITSINTRACK] // input/output
						)
{

 Short_t i,
	itrack,
	ihit,
	j,
	k,
	ntot,
	nadd;

 Double_t angle,
	deltaZ,
	dist,
	dist1,
	Zpos;

 const Double_t NTIMES=0.4;
//	const Double_t NTIMES=1.;


 for(itrack=StartTrackCand; itrack<EndTrackCand; itrack++){
	if( ! keepit[itrack] ) continue;
	if( ! Mvdhits[itrack] ) continue;
	if( Fi_low_limit[itrack] < -99998.) continue;  // case in which
		// the track is completely inside the Mvd region.
	nParHitsinTrack[itrack]=0;

//	loop over the STT parallel hits and try to attach to each candidate track; in this
//	way in one shot I collect also the previously non collected hits and I remove the
//	spurious hits.

	nadd=0;
	for(i=0; i<nSttParHit; i++){
		ihit = fListSttParHits[i];

		if( !SingleHitListStt[ihit] ) continue;
		angle = atan2(info[ihit][1]-fOy[itrack],info[ihit][0]-fOx[itrack]);
		if(angle<0.) angle += 2.*PI;

		// selection on Zpos under assumption that for parallel STT
		// makes 1 turn
		if( fabs(KAPPA[itrack]) > 1.e-20){
			deltaZ = 2.*PI/KAPPA[itrack];
			Zpos = (angle - FI0[itrack])/KAPPA[itrack];

			if(
			    fabs(Zpos-info[ihit][2])>1.5*info[ihit][4] &&
			    fabs(Zpos+deltaZ-info[ihit][2])>1.5*info[ihit][4] &&
			    fabs(Zpos-deltaZ-info[ihit][2])>1.5*info[ihit][4]
			  ) continue;
		}

		if( angle>Fi_up_limit[itrack]){
			angle -= 2.*PI;
			if(angle < Fi_low_limit[itrack]) continue;
		} else if (angle < Fi_low_limit[itrack]){
			angle += 2.*PI;
			if(angle > Fi_up_limit[itrack]) continue;
		}

		dist1 = fabs(
				sqrt(
			(fOx[itrack]-info[ihit][0])*(fOx[itrack]-info[ihit][0])+
			(fOy[itrack]-info[ihit][1])*(fOy[itrack]-info[ihit][1])
				) - fR[itrack]
			);
		dist = fabs( dist1 - info[ihit][3] );


		// info[ihit][3] = drift radius; if  info[ihit][3]>dist then
		// the drift circle of this straw crosses the trajectory.
		if(dist<NTIMES*STRAWRADIUS || info[ihit][3]> dist ){
			ListParHitsinTrack[itrack][ nParHitsinTrack[itrack] ]
				=ihit;
			nParHitsinTrack[itrack]++;
		}


	}	// end of     for(i=0; i<nSttParHit; i++)




 }	// end of for(itrack=StartTrackCand; itrack<EndTrackCand; itrack++)

 return;
}

//------------------------- end of function  PndTrkTracking::CollectParSttHitsagain



//-------------------------  begin of function  PndTrkTracking::EliminateSpuriousSZ

void PndTrkTracking::EliminateSpuriousSZ(
	Short_t MaxTurnofTracks,
	Short_t *nPixelHitsinTrack,
	Short_t *ListPixelHitsinTrack,
	Short_t *nStripHitsinTrack,
	Short_t *ListStripHitsinTrack,
	Short_t *nSkewHitsinTrack,
	Short_t *ListSkewHitsinTrack,
	Double_t *S,
	Double_t *ZED,
	Double_t *DriftRadius,
	Double_t *ErrorDriftRadius,
	Double_t *SchosenPixel,
	Double_t *SchosenStrip,
	Double_t *SchosenSkew,
	Double_t *ZchosenPixel,
	Double_t *ZchosenStrip,
	Double_t *ZchosenSkew,
	Double_t *ErrorchosenPixel,
	Double_t *ErrorchosenStrip,
	Double_t *ErrorchosenSkew,
	Double_t KAPPA,
	Double_t FI0,
	Double_t Rr
	)
{

	bool	already[MAXSTTHITS];

	Short_t i,
		 j,
		 k,
		 Kmin,
		 auxnMvdPixel,
		 auxListMvdPixel[MAXMVDPIXELHITS],
		 auxnMvdStrip,
		 auxListMvdStrip[MAXMVDSTRIPHITS],
		 auxnSttSkew,
		 auxListSttSkew[MAXSTTHITS];

	Int_t	len,
		 nr2,
		 nrounds0,
		 nrounds1,
		 nchosen,
		 Nround[4];

	const Double_t  MvdCut=0.8,
//	const Double_t  MvdCut=0.3,
//	const Double_t  MvdCut=0.1,
			minimumSttDriftError = 1.;

	Double_t ddd,
		 dista,
		 dista1,
		 dista0,
		 d_min,
		 error,
		 zeta0,
		 zeta1,
		 Dista[4],
		 dista_storage[MAXSTTHITS],
		 Errore[4],
		 Esse[4],
		 Zeta[4];

 PndTrkCTGeometryCalculations GeomC;

	len = sizeof(already);
	memset (already,false,len);

	auxnMvdPixel=0;
	auxnMvdStrip=0;
	auxnSttSkew=0;


	for(i=0;i<*nPixelHitsinTrack;i++){
		if( fabs(GeomC.Dist_SZ(Rr,KAPPA,FI0,ZED[i],S[i],&nrounds0)) < MvdCut
						&&
					abs(nrounds0)<=MaxTurnofTracks){
			auxListMvdPixel[auxnMvdPixel]=ListPixelHitsinTrack[i];
			if(nrounds0>=0)SchosenPixel[ListPixelHitsinTrack[i]]=S[i]+nrounds0*2.*PI;
			else SchosenPixel[ListPixelHitsinTrack[i]]=S[i]+(nrounds0-1)*2.*PI;
			ZchosenPixel[ListPixelHitsinTrack[i]]=ZED[i];
			ErrorchosenPixel[ListPixelHitsinTrack[i]]=ERRORPIXEL;
			auxnMvdPixel++;
		}
	}	// end of  for(i=0;i<*nPixelHitsinTrack;i++)

	for(j=0;j<*nStripHitsinTrack;j++){
		i=j+(*nPixelHitsinTrack);

		if( fabs(GeomC.Dist_SZ(Rr,KAPPA,FI0,ZED[i],S[i],&nrounds0)) < MvdCut
						&&
					abs(nrounds0)<=MaxTurnofTracks){
			auxListMvdStrip[auxnMvdStrip]=ListStripHitsinTrack[j];
			if(nrounds0>=0)SchosenStrip[ListStripHitsinTrack[j]]=S[i]+nrounds0*2.*PI;
			else SchosenStrip[ListStripHitsinTrack[j]]=S[i]+(nrounds0-1)*2.*PI;
			ZchosenStrip[ListStripHitsinTrack[j]]=ZED[i];
			ErrorchosenStrip[ListStripHitsinTrack[j]]=ERRORSTRIP;
			auxnMvdStrip++;
		}
	}	// end of  for(j=0;j<*nStripHitsinTrack;j++)

//-----------------stampe.
if(istampa>3){
	cout<<"in EliminateSpuriousSZ : nSkew hit = "<<*nSkewHitsinTrack
	<<", K = "<< KAPPA<<", FI0 = "<< FI0 <<endl;
}
//--------------fine stampe.

	for(j=0;j<*nSkewHitsinTrack;j++){
		i=j+(*nPixelHitsinTrack)+
		    (*nStripHitsinTrack);


//-----------------stampe.
if(istampa>3){
	cout<<"in EliminateSpuriousSZ : skew hit n. "<<ListSkewHitsinTrack[j];
	cout<<"  ZED[i] = "<<ZED[i]<<", S[i] "<<S[i];
	cout<<endl;
}
//--------------fine stampe.
			dista =
		fabs(GeomC.Dist_SZ(Rr,KAPPA,FI0,ZED[i]+DriftRadius[i],S[i],&nrounds0));
			ddd = fabs(GeomC.Dist_SZ(Rr,KAPPA,FI0,ZED[i]-DriftRadius[i],S[i],&nr2));
			if( abs(nrounds0) > MaxTurnofTracks &&
				 abs(nr2)>MaxTurnofTracks)
			{
				continue;
			} else if ( abs(nr2)>MaxTurnofTracks)
			{
				 zeta0 = ZED[i]+DriftRadius[i];
				 nchosen=nrounds0;
			} else if ( abs(nrounds0) > MaxTurnofTracks )
			{
				 dista = ddd;
				 zeta0 = ZED[i]-DriftRadius[i];
				 nchosen=nr2;
			} else
			{
				if( fabs(dista)>fabs(ddd) ) {
				 dista = ddd;
				 zeta0 = ZED[i]-DriftRadius[i];
				 nchosen=nr2;
				} else {
				 zeta0 = ZED[i]+DriftRadius[i];
				 nchosen=nrounds0;
				}
			}
			error = DriftRadius[i];
			if(nchosen>=0)SchosenSkew[ListSkewHitsinTrack[j]]=S[i]+nchosen*2.*PI;
			else SchosenSkew[ListSkewHitsinTrack[j]]=S[i]+(nchosen-1)*2.*PI;
			ZchosenSkew[ListSkewHitsinTrack[j]]=zeta0;
//-----------------stampe.
if(istampa>3){
	cout<<"in EliminateSpuriousSZ : insomma, dista prima della selezione = "<<
	dista<< ", ed e' da comparare\n\tcon 4*error = "<<  4.*error<<endl<<
	"\toppure con 2.*minimumSttDriftError = "<<2.*minimumSttDriftError<<endl;
}
//--------------fine stampe.


		if(	dista < 4.*error
//			dista < 1.1*error
				||
			dista < 2.*minimumSttDriftError
			)
		{

			//  check now if the other solution of the same hit
			//  has already been selected before;

			if( already[ ListSkewHitsinTrack[j] ]){
			  // in this case the other solution has already been
			  // selected; therefore in   dista_storage[ ListSkewHitsinTrack[j] ]
			  // there must be the distance  previously calculated;
			  if( dista_storage[ ListSkewHitsinTrack[j] ]>dista ){
				auxListSttSkew[auxnSttSkew]=ListSkewHitsinTrack[j];
				ErrorchosenSkew[ListSkewHitsinTrack[j]]=error;
				auxnSttSkew++;
			  }
			} else {  // continuation of  if( already[ ListSkewHitsinTrack[j] ] 


			  already[ ListSkewHitsinTrack[j] ] = true;
			  // store the  distance calculated;
			  dista_storage[ListSkewHitsinTrack[j] ]=dista;
			  auxListSttSkew[auxnSttSkew]=ListSkewHitsinTrack[j];
			  ErrorchosenSkew[ListSkewHitsinTrack[j]]=error;
			  auxnSttSkew++;




			} // end of   if( already[ ListSkewHitsinTrack[j] ] )





		}  // end of  if( dista < 4.*error ....)
	}	// end of  for(j=0;j<*nSkewHitsinTrack;j++)




//	reload the list of good hits

	*nPixelHitsinTrack = auxnMvdPixel;
	*nStripHitsinTrack = auxnMvdStrip;
	*nSkewHitsinTrack = auxnSttSkew;
	for(j=0;j<*nPixelHitsinTrack;j++){
		ListPixelHitsinTrack[j]=auxListMvdPixel[j];
	}	// end of  for(j=0;j<*nPixelHitsinTrack;j++)
	for(j=0;j<*nStripHitsinTrack;j++){
		ListStripHitsinTrack[j]=auxListMvdStrip[j];
	}	// end of  for(j=0;j<*nStripHitsinTrack;j++)
	for(j=0;j<*nSkewHitsinTrack;j++){
		ListSkewHitsinTrack[j]=auxListSttSkew[j];
	}	// end of  for(j=0;j<*nSkewHitsinTrack;j++)


	return;

}
//-------------------------  end of function  PndTrkTracking::EliminateSpuriousSZ




//----------------  begin of function PndTrkTracking::ExtractInfoFromMvdTrackCand

void PndTrkTracking::ExtractInfoFromMvdTrackCand()
{

  Short_t
	i,
	j,
	k;

  Double_t qop;

  PndTrackCand* pMvdTrackCand;

  PndTrackCandHit pndtrackcandhit; 

  for( i= 0; i< fnMvdTrackCand ; i++){

	pMvdTrackCand= (PndTrackCand*) fMvdTrackCandArray->At(i);
	TVector3 dirSeed=pMvdTrackCand->getDirSeed();
	TVector3 posSeed=pMvdTrackCand->getPosSeed();
	qop = pMvdTrackCand->getQoverPseed();
	// n. hits in this track cand
	fnHitMvdTrackCand[i] = pMvdTrackCand->GetNHits();
	if( fnHitMvdTrackCand[i]>MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK){
		cout<<"from PndTrkTracking, fnHitMvdTrackCand[i] = "<<fnHitMvdTrackCand[i]
		<<" and it is > MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK (="
		<<MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK
		<<"); setting fnHitMvdTrackCand[i] to "
		<<MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK<<endl;
	fnHitMvdTrackCand[i]=MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK;
	}

	Short_t kPixel,kStrip;

	for(j=0, k=0, kPixel=0, kStrip=0; j<fnHitMvdTrackCand[i]; j++){
		pndtrackcandhit = pMvdTrackCand->GetSortedHit(j);

	// the following case should never happen (in principle), but, just to be on
	// the safe side ....
		if(pndtrackcandhit.GetHitId()<0||pndtrackcandhit.GetDetId()<0)continue;

	// this is a Pixel.
	if( pndtrackcandhit.GetDetId()==
		FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)
			&& kPixel < MAXMVDPIXELHITSINTRACK){

		// the following is a protection because the maximum Mvd hit n.
		// cannot exceed MAXMVDPIXELHITS .
		if( pndtrackcandhit.GetHitId()>MAXMVDPIXELHITS ){
			cout<<"from PndTrkTracking, this Pixel Mvd hit has a number = "
			<<pndtrackcandhit.GetHitId()<<
			" that is > MAXMVDPIXELHITS (="<<
				MAXMVDPIXELHITS<<"), rejected!"<<endl;
			continue;
		}


		finMvdTrackCandPixel[ pndtrackcandhit.GetHitId() ]= true;
		fListHitTypeMvdTrackCand[i][k] = pndtrackcandhit.GetDetId();
					  // this in reality is the Branch name of
					  // Hit; this Branch name is used to identify
					  // a Pixel.
					  // If it is -1 is supposed to be noise but
					  // a in which Pixel or Strip?? Mistery.
		kPixel++;

	// this is a Strip.
	} else if( pndtrackcandhit.GetDetId()==
		FairRootManager::Instance()->GetBranchId(fMvdStripBranch)
			&& kStrip < MAXMVDSTRIPHITSINTRACK){

		// the following is a protection because the maximum Mvd hit n. cannot
		// exceed MAXMVDSTRIPHITS .
		if( pndtrackcandhit.GetHitId()>MAXMVDSTRIPHITS ){
			cout<<"from PndTrkTracking, this Strip Mvd hit has a number = "
			<<pndtrackcandhit.GetHitId()<<
			" that is > MAXMVDSTRIPHITS (="<<
			MAXMVDSTRIPHITS<<"), rejected!"<<endl;
			continue;
		}

		finMvdTrackCandStrip[ pndtrackcandhit.GetHitId() ]= true;
		fListHitTypeMvdTrackCand[i][k] = pndtrackcandhit.GetDetId();
					  // this in reality is the Branch name of
					  // Hit; this Branch name is used to identify
					  // a Pixel.
					  // If it is -1 is supposed to be noise but
					  // a in which Pixel or Strip?? Mistery.
		kStrip++;
	} else {	// this is the case should (in principle) never happen.
		continue;	// ignore this hit.
	}

	// the following is the native Hit number that one can use
	// to extract all the info. If it is -1 I won't consider it (noise).
	fListHitMvdTrackCand[i][k] = pndtrackcandhit.GetHitId();

	k++;

     }    //   end of    for(j=0; j<fnHitMvdTrackCand; j++)

     fnHitMvdTrackCand[i]=k; // if the case, readjust the # of hits in this candidate.

  }   //   end of        for( i= 0; i< fnMvdTrackCand ; i++)


//  now load the arrays containing the Pixel Hits and Strip Hits
// not belonging to any Mvd Track Cand

  fnMvdUSPixelHitNotTrackCand=0;
  fnMvdDSPixelHitNotTrackCand=0;
  fnMvdUSStripHitNotTrackCand=0;
  fnMvdDSStripHitNotTrackCand=0;
  for( i= 0; i< fnMvdPixelHit ; i++){
	if( ! finMvdTrackCandPixel[i] ){
		if( fZMvdPixel[i]>=0.){
			fListMvdDSPixelHitNotTrackCand[fnMvdDSPixelHitNotTrackCand] = i;
			fnMvdDSPixelHitNotTrackCand++;
		} else {
			fListMvdUSPixelHitNotTrackCand[fnMvdUSPixelHitNotTrackCand] = i;
			fnMvdUSPixelHitNotTrackCand++;
		}
	}
  }

  for( i= 0; i< fnMvdStripHit ; i++){
	if( ! finMvdTrackCandStrip[i] ){
		if( fZMvdStrip[i]>=0.){
			fListMvdDSStripHitNotTrackCand[fnMvdDSStripHitNotTrackCand] = i;
			fnMvdDSStripHitNotTrackCand++;
		} else {
			fListMvdUSStripHitNotTrackCand[fnMvdUSStripHitNotTrackCand] = i;
			fnMvdUSStripHitNotTrackCand++;
		}
	}
  }

 
 }
//-----------------------  end of function   PndTrkTracking::ExtractInfoFromMvdTrackCand


//----------begin of function PndTrkTracking::FindCharge

void PndTrkTracking::FindCharge(
	Double_t oX,
	Double_t oY,
	Short_t nParallelHits,
	Double_t *X,
	Double_t *Y,
	Short_t  * Charge
	)
{

	Short_t ihit,
		nleft,
		nright;

	Double_t cross,
		 disq,
		 minl,
		 minr;


	// this methods works with the hypothesis that this track comes
	//  from (0,0)

	for(ihit=0, nleft=0, nright=0, minr = 9999999., minl = 9999999.; ihit<nParallelHits; ihit++){ 
	// find the Z component of the cross product between the vector from (0,0) to center of
	// circular trajectory [namely, (oX,oY) ]  and the Position vector of the center of the
	// parallel Hits [namely, (x,y)].

		cross = oX*Y[ihit] -
			oY*X[ihit];

	// if  cross >0  hits stays 'on the left' (which means clockwise to go from the origin
	// to the hit following the smaller path) otherwise it stays 'on the right'.

		if (cross>0.) {
			disq =	X[ihit]*X[ihit]+Y[ihit]*Y[ihit];
			nleft++;
		} else {
			nright++;
		}
	}	// end of   for(ihit=0, nleft=0, nright=0;....

	if( nright> nleft) {
		*Charge = -1;
	} else if ( nleft > nright) {
		*Charge = 1;
	} else {	// then choose according the closest hit ti the center
		if( minr < minl ) *Charge = -1;
		else  *Charge = 1;
	}



}

//----------end of function PndTrkTracking::FindCharge





//----------start  function PndTrkTracking::FindingParallelTrackAngularRange

void  PndTrkTracking::FindingParallelTrackAngularRange(
	Double_t oX,
	Double_t oY,
	Double_t Rr,
	Short_t  Charge,
	Double_t *Fi_low_limit,
	Double_t *Fi_up_limit,
	Short_t * status,
	Double_t Rmi,	// Rmin of cylindrical volume intersected by track;
	Double_t Rma	// Rmax of cylindrical volume intersected by track;
	)
{
// -------------- calculate the maximum fi and minimum fi spanned by this track,

// see logbook pag.270; by using the Rmin and Rmax of the straw detector.
// this function works also when circular trajectory in XY doesn't pass through (0,0).

	bool	intersection_inner,
		intersection_outer;
	Double_t	teta1,
			teta2,
			tetavertex,
			a,
			cosT,
			cost,
			cosFi,
			cosfi,
			Fi,
			fi,
			FI0,
			Px,
			Py,
			tmp;


	Rma += 1. ; // add a safety margin.
	Rmi -= 1. ; // add a safety margin.




	a = sqrt(oX*oX+oY*oY);

	//  preliminary condition
	if(a + Rr <= Rmi )	// in this case there might be hits at radius < Rmi.
		 { *status = -1 ;return;}
	if( a >= Rr + Rma || Rr >= a + Rma)  // in this case there can be no hits with radius < Rma.
		 { *status = -2;return;}

	if( a - Rr >= Rmi ) intersection_inner = false; else intersection_inner = true;

	if( a + Rr <= Rma || a - Rr >= Rma  )
		 intersection_outer = false; else intersection_outer = true;

	if( (! intersection_inner) && (! intersection_outer) ){
		*Fi_low_limit = 0.;
		*Fi_up_limit = 2.*PI;
		*status = 1;
		return;
	}

//	now the calculation

	FI0 = atan2(-oY,-oX);
	if( intersection_outer ){
		cosFi = (a*a + Rr*Rr - Rma*Rma)/(2.*Rr*a);
		if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.;
		Fi = acos(cosFi);
	}

	if( intersection_inner ){
		cosfi = (a*a + Rr*Rr - Rmi*Rmi)/(2.*Rr*a);
		if(cosfi<-1.) cosfi=-1.; else if(cosfi>1.) cosfi=1.;
		fi = acos(cosfi);
	}


	if( Charge < 0.){ // this particle rotates counterclockwise when looking into the beam
		if( intersection_outer && intersection_inner){
			*Fi_low_limit=FI0 + fi;
			*Fi_up_limit= FI0 +Fi;

		} else if (intersection_inner) {
			*Fi_low_limit=FI0 + fi;
			*Fi_up_limit= FI0 - fi;
		} else {
			*Fi_low_limit=FI0 - Fi;
			*Fi_up_limit= FI0 + Fi;
		}	// end of    if( intersection_outer && intersection_inner



	} else {	// continuation of   if( Charge < 0.)

		if( intersection_outer && intersection_inner){
			*Fi_low_limit=FI0 - Fi;
			*Fi_up_limit= FI0 - fi;

		} else if (intersection_inner) {
			*Fi_low_limit=FI0 + fi;	// must invert because low limit must be < up limit
			*Fi_up_limit= FI0 - fi;
		} else {
			*Fi_low_limit=FI0 - Fi;
			*Fi_up_limit= FI0 + Fi;
		}	// end of    if( intersection_outer && intersection_inner


	}	// end of  if( Charge < 0.)




	if(*Fi_low_limit<0.) {
		*Fi_low_limit=fmod(*Fi_low_limit,2.*PI);
		*Fi_low_limit += 2.*PI;
	} else if (*Fi_low_limit>=2.*PI){
		*Fi_low_limit=fmod(*Fi_low_limit,2.*PI);
	}
	if(*Fi_up_limit<0.) {
		*Fi_up_limit=fmod(*Fi_up_limit,2.*PI);
		*Fi_up_limit += 2.*PI;
	} else if (*Fi_up_limit>=2.*PI){
		*Fi_up_limit=fmod(*Fi_up_limit,2.*PI);
	}

	//	Modify *Fi_up_limit by adding
	//	2PI if it is the case, in order to make *Fi_up_limit > *Fi_low_limit.
	if( *Fi_up_limit < *Fi_low_limit ) *Fi_up_limit += 2.*PI;
	if( *Fi_up_limit < *Fi_low_limit ) *Fi_up_limit = *Fi_low_limit;

	*status = 0;



      return;
}


//---------- end of  function PndTrkTracking::FindingParallelTrackAngularRange



//----------start of function PndTrkTracking::FixDiscontinuitiesFiangleinSZplane

void PndTrkTracking::FixDiscontinuitiesFiangleinSZplane(
	Short_t TemporarynSkewHitsinTrack,
	Vec <Double_t> & S,
	Double_t *Fi_initial_helix_referenceframe,
	Short_t Charge
	)
{

 Short_t i;
 Double_t max, min;

 if( Charge >0 )
 {
	for(i=0 ; i<TemporarynSkewHitsinTrack; i++){
		if( S[i] > *Fi_initial_helix_referenceframe )  S[i]-= 2.*PI;
	}
 } else {
	for(i=0 ; i<TemporarynSkewHitsinTrack; i++){
		if( S[i] < *Fi_initial_helix_referenceframe )  S[i]+= 2.*PI;
	}
 }

     return;

}
//----------end of function PndTrkTracking::FixDiscontinuitiesFiangleinSZplane


//---------- begin of function PndTrkTracking::InfoXYZParal

void PndTrkTracking::InfoXYZParal(
	Double_t info[][7],
	Short_t infopar,
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t KAPPA,
	Double_t FI0,
	Short_t Charge,
	Double_t *Posiz
	)
{

   Double_t fi, norm, vers[2];

   vers[0] = Oxx - info[infopar][0];
   vers[1] = Oyy - info[infopar][1];
   norm = sqrt( vers[0]*vers[0] + vers[1]*vers[1] );

   if(norm < 1.e-20) {
     Posiz[0] = -999999999.;
     return;
   }
   





   if( fabs( Rr - fabs( norm - info[infopar][3] ) ) // distance trajectory-drift radius
				<
		fabs( Rr - (norm + info[infopar][3]) )  ) {

	Posiz[0] = info[infopar][0] + info[infopar][3]*vers[0]/norm;
	Posiz[1] = info[infopar][1] + info[infopar][3]*vers[1]/norm;

   } else {

	Posiz[0] = info[infopar][0] - info[infopar][3]*vers[0]/norm;
	Posiz[1] = info[infopar][1] - info[infopar][3]*vers[1]/norm;

   }	// end of if ( fabs( Rr - fabs( Distance - info[infopar][3] ) ).....




//   Posiz[0] = info[infopar][0] + info[infopar][3]*vers[0]/norm;
//   Posiz[1] = info[infopar][1] + info[infopar][3]*vers[1]/norm;

   if( fabs(KAPPA)<1.e-20 ){
     Posiz[2] = -888888888.;
     return;
   }


   fi = atan2(-vers[1],-vers[0]);
   if(fi<0.)  fi += 2.*PI;

   if ( Charge > 0){
    if(fi > FI0 )  FI0 += 2.*PI;
//    Posiz[2] = (FI0-fi)/KAPPA;
   } else {
    if(fi < FI0 )  fi += 2.*PI;
   }
    Posiz[2] = (fi-FI0)/KAPPA;

   return;
}

//----------end of function PndTrkTracking::InfoXYZParal



//--------begin of function PndTrkTracking::Initial_SttParHits_DecreasingR_Ordering

void  PndTrkTracking::Initial_SttParHits_DecreasingR_Ordering(
	Double_t info[][7],
	Short_t *ListSttParHi,
	Int_t nSttParHit
	)
{

 Short_t
	j,
	OLDListSttParHits[nSttParHit];

 Int_t
	auxIndex[nSttParHit];

 Double_t
	auxRvalues[nSttParHit];


//   ordering the parallel hits by decreasing spatial radius.

 for (j = 0; j< nSttParHit; j++){
	auxIndex[j]=j;
	auxRvalues[j]=
		info[ListSttParHi[ j ]  ][0]*info[ListSttParHi[ j ]  ][0]+
		info[ListSttParHi[ j ]  ][1]*info[ListSttParHi[ j ]  ][1];
	OLDListSttParHits[j]=ListSttParHi[j];
 }

 PndTrkMergeSort Sorter;

 Sorter.Merge_Sort( (Short_t) nSttParHit, auxRvalues, auxIndex);

 for (j = 0; j< nSttParHit; j++){
	ListSttParHi[ nSttParHit-1-j] = OLDListSttParHits[ auxIndex[ j ]   ];
 }

 return;
}

//--------end function PndTrkTracking::Initial_SttParHits_DecreasingR_Ordering




//----------------------  begin function   PndTrkTracking::LoadPndTrack_TrackCand

void PndTrkTracking::LoadPndTrack_TrackCand(
	Vec <bool>& keepit,
	Vec <bool>& SttSZfit,
	Short_t nTotalCandidates,
	Vec <Short_t>& Charge,
	Int_t nSttTrackCand,
	Double_t *FI0,
	Double_t *KAPPA,
	Double_t info[][7],
	Double_t SchosenSkew[][MAXSTTHITS],
	Double_t ZchosenSkew[][MAXSTTHITS]
			)
{

//-------  load the new PndTrackCand ; each track has the STT and the Mvd hits associated
//-------  also load the new PndTrack ; each track has the STT and the Mvd hits associated

 Int_t
	ipinco,
	j,
	k,
	l,
	ncand
	;

 Double_t
	ddd,
	dis,
	Distance,
	Oxx,
	Oyy,
	Ptras,
	Pxini,
	Pyini,
	Pzini,
	px,
	py,
	qop,
	x,
	y,
	Posiz1[3],
	versor[2]
	;

 TVector3
	ErrMomentum,
	ErrPosition,
	Momentum,
	Position
	;

//--------------------------
int nSttHit = fSttHitArray->GetEntriesFast();

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




 for(ncand=0, ipinco = 0; ncand< nTotalCandidates; ncand++){
	if(!keepit[ncand]) continue;
	// case in which there was no Skew hits and no KAPPA info and that
	// candidate could not be associated to any Mvd hits --> no KAPPA information!
	if(ncand<nSttTrackCand && ! SttSZfit[ncand]) continue;
	Oxx = fOx[ncand];
	Oyy = fOy[ncand];
	dis=sqrt( Oxx*Oxx+Oyy*Oyy );
	if( dis < 1.e-20)  continue;
	Ptras = fR[ncand]*0.003*BFIELD;
	Pxini = -Charge[ncand]*Ptras*Oyy/dis;
	Pyini = Charge[ncand]*Ptras*Oxx/dis;

//   starting point not necessarily at x=0., y=0.

	x= fOx[ncand] + fR[ncand]*cos(FI0[ncand]);
	y= fOy[ncand] + fR[ncand]*sin(FI0[ncand]);
	TVector3 posSeed(x,y,0.);  //  the starting point

	if(fabs(KAPPA[ncand])>1.e-20  ){
		Pzini = -Charge[ncand]*0.003*BFIELD/KAPPA[ncand];
		if(fabs(Pzini) > PMAX)  continue;
	} else {
		continue;
	}
	// PndTrackCand Array loading
	new((*fSttMvdPndTrackCandArray)[ipinco])  PndTrackCand;
	PndTrackCand *pTrckCand = (PndTrackCand*) fSttMvdPndTrackCandArray->At(ipinco);
	TVector3 dirSeed(Pxini,Pyini,Pzini); // momentum direction in starting point
	qop = Charge[ncand]/dirSeed.Mag();
	dirSeed.SetMag(1.);
	pTrckCand->setTrackSeed(posSeed, dirSeed, qop);
	pTrckCand->setMcTrackId(  -1   );

	for(j=0; j< fnTrackCandHit[ncand]; j++){

	     switch (fListTrackCandHitType[ncand][j]){
		case 0:
			pTrckCand->AddHit(FairRootManager::Instance()->
				GetBranchId(fMvdPixelBranch),
				(Int_t)fListTrackCandHit[ncand][j],j);
		break;
		case 1:
			pTrckCand->AddHit(FairRootManager::Instance()->
				GetBranchId(fMvdStripBranch),
				(Int_t)fListTrackCandHit[ncand][j],j);
		break;
		case 2:
			pTrckCand->AddHit(FairRootManager::Instance()->
				GetBranchId(fSttBranch),
				(Int_t)fListTrackCandHit[ncand][j],j);
		break;
		case 3:
			pTrckCand->AddHit(FairRootManager::Instance()->
				GetBranchId(fSttBranch),
				(Int_t)fListTrackCandHit[ncand][j],j);
		break;
	     }
	}

	// PndTrack Array loading
	//  last = last hit in found track; first = first hit in found track;
	//  pTrckCand =  pointer to the corresponding  PndTrackCand.

	//  the first hit

	if (fListTrackCandHitType[ncand][0] == 0){  //  Mvd Pixel
		Posiz1[0] = fXMvdPixel[ fListTrackCandHit[ncand][0] ];
		Posiz1[1] = fYMvdPixel[ fListTrackCandHit[ncand][0] ];
		Posiz1[2] = fZMvdPixel[ fListTrackCandHit[ncand][0] ];
		ErrPosition.SetX(fsigmaXMvdPixel[ fListTrackCandHit[ncand][0] ]/sqrt(12.));
		ErrPosition.SetY(fsigmaXMvdPixel[ fListTrackCandHit[ncand][0] ]/sqrt(12.));
		ErrPosition.SetZ(fsigmaXMvdPixel[ fListTrackCandHit[ncand][0] ]/sqrt(12.));
	} else if (fListTrackCandHitType[ncand][0] == 1){  //  Mvd Strip
		Posiz1[0] = fXMvdStrip[ fListTrackCandHit[ncand][0] ];
		Posiz1[1] = fYMvdStrip[ fListTrackCandHit[ncand][0] ];
		Posiz1[2] = fZMvdStrip[ fListTrackCandHit[ncand][0] ];
		ErrPosition.SetX(fsigmaXMvdStrip[ fListTrackCandHit[ncand][0] ]/sqrt(12.));
		ErrPosition.SetY(fsigmaXMvdStrip[ fListTrackCandHit[ncand][0] ]/sqrt(12.));
		ErrPosition.SetZ(fsigmaXMvdStrip[ fListTrackCandHit[ncand][0] ]/sqrt(12.));
	} else if( fListTrackCandHitType[ncand][0] == 2 ){  // it is a parallel straw hit
		InfoXYZParal(
			info,
			fListTrackCandHit[ncand][0],
			fOx[ncand],
			fOy[ncand],
			fR[ncand],
			KAPPA[ncand],
			FI0[ncand],
			Charge[ncand],
			Posiz1
			);

		// cases in which the calculation of the position failed, see
		// InfoXYZParal  method.
		if(Posiz1[2]<-888888887. || Posiz1[0] < -999999998.)
			continue;

		ErrPosition.SetX(0.02);	// 200 microns
		ErrPosition.SetY(0.02);	// 200 microns
		ErrPosition.SetZ(1.);		// 1 cm


	} else if ( fListTrackCandHitType[ncand][0] == 3 ){  // it is a skew straw hit

		Posiz1[0] = fOx[ncand]+fR[ncand]*cos(SchosenSkew[ncand][ fListTrackCandHit[ncand][0] ]);
		Posiz1[1] = fOy[ncand]+fR[ncand]*sin(SchosenSkew[ncand][ fListTrackCandHit[ncand][0] ]);
		Posiz1[2] = ZchosenSkew[ncand][ fListTrackCandHit[ncand][0] ];
		ErrPosition.SetX(0.02);	// 200 microns
		ErrPosition.SetY(0.02);	// 200 microns
		ErrPosition.SetZ(1.);		// 1 cm

	}

	Position.SetX( Posiz1[0] );
	Position.SetY( Posiz1[1] );
	Position.SetZ( Posiz1[2] );
	versor[0] = fOx[ncand]-Posiz1[0];
	versor[1] = fOy[ncand]-Posiz1[1];
	Distance = sqrt(versor[0]*versor[0]+versor[1]*versor[1]);
	versor[0] /= Distance;
	versor[1] /= Distance;
	px = -Charge[ncand]*Ptras*versor[1];
	py = Charge[ncand]*Ptras*versor[0];
	Momentum.SetX(px);
	Momentum.SetY(py);
	Momentum.SetZ(Pzini);
	ErrMomentum.SetX(0.05*Ptras); //  set at 5% all the times.
	ErrMomentum.SetY(0.05*Ptras); //  set at 5% all the times.
	ErrMomentum.SetZ(0.05*Pzini); //  set at 5% all the times.
	//  the plane of this FairTrackParP better is perpendicular to
	//  the momentum direction
	ddd = Ptras*sqrt(Ptras*Ptras+Pzini*Pzini);

	FairTrackParP first( Position,  Momentum,
	   ErrPosition, ErrMomentum, Charge[ncand],
	   Position,
	   TVector3(py/Ptras, -px/Ptras, 0.), //  first vector defining the plane
	   TVector3(Pzini*px/ddd,Pzini*py/ddd,-Ptras*Ptras/ddd) //second vector defining the plane
				);
	//  the last hit

	k = fnTrackCandHit[ncand]-1;
	if (fListTrackCandHitType[ncand][k] == 0){  //  Mvd Pixel
		Posiz1[0] = fXMvdPixel[ fListTrackCandHit[ncand][k] ];
		Posiz1[1] = fYMvdPixel[ fListTrackCandHit[ncand][k] ];
		Posiz1[2] = fZMvdPixel[ fListTrackCandHit[ncand][k] ];
		ErrPosition.SetX(fsigmaXMvdPixel[ fListTrackCandHit[ncand][k] ]/sqrt(12.));
		ErrPosition.SetY(fsigmaXMvdPixel[ fListTrackCandHit[ncand][k] ]/sqrt(12.));
		ErrPosition.SetZ(fsigmaXMvdPixel[ fListTrackCandHit[ncand][k] ]/sqrt(12.));
	} else if (fListTrackCandHitType[ncand][k] == 1){  //  Mvd Strip
		Posiz1[0] = fXMvdStrip[ fListTrackCandHit[ncand][k] ];
		Posiz1[1] = fYMvdStrip[ fListTrackCandHit[ncand][k] ];
		Posiz1[2] = fZMvdStrip[ fListTrackCandHit[ncand][k] ];
		ErrPosition.SetX(fsigmaXMvdStrip[ fListTrackCandHit[ncand][k] ]/sqrt(12.));
		ErrPosition.SetY(fsigmaXMvdStrip[ fListTrackCandHit[ncand][k] ]/sqrt(12.));
		ErrPosition.SetZ(fsigmaXMvdStrip[ fListTrackCandHit[ncand][k] ]/sqrt(12.));
	} else if( fListTrackCandHitType[ncand][k] == 2 ){  // it is a parallel straw hit
		InfoXYZParal (
			info,
			fListTrackCandHit[ncand][k],
			fOx[ncand],
			fOy[ncand],
			fR[ncand],
			KAPPA[ncand],
			FI0[ncand],
			Charge[ncand],
			Posiz1
			);
		if(Posiz1[2]<-888888887. || Posiz1[0] < -999999998.)
			continue;

		ErrPosition.SetX(0.02);	// 200 microns
		ErrPosition.SetY(0.02);	// 200 microns
		ErrPosition.SetZ(1.);	// 1 cm

	} else if ( fListTrackCandHitType[ncand][k] == 3 ){  // it is a skew straw hit

		Posiz1[0] = fOx[ncand]+fR[ncand]*cos(SchosenSkew[ncand][ fListTrackCandHit[ncand][k] ]);
		Posiz1[1] = fOy[ncand]+fR[ncand]*sin(SchosenSkew[ncand][ fListTrackCandHit[ncand][k] ]);
		Posiz1[2] = ZchosenSkew[ncand][ fListTrackCandHit[ncand][k] ];
		ErrPosition.SetX(0.02);	// 200 microns
		ErrPosition.SetY(0.02);	// 200 microns
		ErrPosition.SetZ(1.);	// 1 cm

	}

	Position.SetX( Posiz1[0] );
	Position.SetY( Posiz1[1] );
	Position.SetZ( Posiz1[2] );
	versor[0] = fOx[ncand]-Posiz1[0];
	versor[1] = fOy[ncand]-Posiz1[1];
	Distance = sqrt(versor[0]*versor[0]+versor[1]*versor[1]);
	versor[0] /= Distance;
	versor[1] /= Distance;
	px = -Charge[ncand]*Ptras*versor[1];
	py = Charge[ncand]*Ptras*versor[0];
	Momentum.SetX(px);
	Momentum.SetY(py);
	Momentum.SetZ(Pzini);
	//  the plane of this FairTrackParP better is perpendicular to
	//  the momentum direction
	ddd = Ptras*sqrt(Ptras*Ptras+Pzini*Pzini);



	FairTrackParP last( Position,  Momentum,
	   ErrPosition, ErrMomentum, Charge[ncand],
	   Position,
	   TVector3(py/Ptras, -px/Ptras, 0.), //  first vector defining the plane
	   TVector3(Pzini*px/ddd,Pzini*py/ddd,-Ptras*Ptras/ddd) //second vector defining the plane
				);



	// loading actually the PndTrack
	PndTrack *pTrck = new((*fSttMvdPndTrackArray)[ipinco]) PndTrack(first,last,*pTrckCand);
	pTrck->SetRefIndex(ipinco);
	pTrck->SetFlag(0);

	ipinco++;
 }	// end of     for(ncand=0, ipinco = 0; ncand< nTotalCandidates; ncand++)


}


//----------end function PndTrkTracking::LoadPndTrack_TrackCand



//----------begin function PndTrkTracking::LoadSZetc_forSZfit


void PndTrkTracking::LoadSZetc_forSZfit(
	Short_t ncand,	// input
	Short_t nhitsinfit,
	Double_t * TemporaryS,		// input
	Double_t * TemporaryZ,		// input
	Double_t * TemporaryZDrift,	// input
	Double_t * TemporaryZErrorafterTilt,	// input
	bool YesGLPKfitSZ,		// input

	Vec <Double_t>& ErrorDriftRadius,	 // output
	Double_t * ErrorDriftRadiusbis,	 // output
	Vec <Double_t>& DriftRadius,		 // output
	Double_t * DriftRadiusbis,	 // output
	Vec <Double_t> & S,			 // output
	Double_t * Sbis,		 // output
	Vec <Double_t>& ZED,			 // output
	Double_t * ZEDbis		 // output
	)
{
 //---------------------   here calculate the S and Z values of Mvd Pixels, Mvd Strips,
 //	 Stt Skew hits and SciTil hits (if present).

 //  the difference between S and Sbis, ZED and ZEDbis, DriftRadius and DriftRadiusbis,
 //  ErrorDriftRadius and ErrorDriftRadiusbis, is that S, ZED etc. contain the list of
 //  Pixel+Strips+SciTil + other Skew Stt hits in case Pixel+Strips+SciTil are <= 2; instead
 //  Sbis, ZEDbis etc. contain Pixel+Strips+all Skew Stt hits.


	int
	i,
	j,
	k,
	kall,
	location;


	// the Mvd Pixels hit
	for(i=0; i< fnMvdPixelHitsinTrack[ncand]; i++){
		k=fListMvdPixelHitsinTrack[ncand][i];
		ZEDbis[i] = ZED[i] = fZMvdPixel[k];
		S[i] = atan2( fYMvdPixel[k]-fOy[ncand],fXMvdPixel[k]-fOx[ncand]);
		if(S[i]<0.) S[i] +=2.*PI;
		Sbis[i] = S[i];
		// DriftRadius is set conventionally at -1, for later use in the SZ fit;
		// the error on the point used in the fit is ErrorDriftRadius and this
		// is overestimated to be  1cm.
		DriftRadiusbis[i]=DriftRadius[i]=-1.;
		if(YesGLPKfitSZ){
		// the following is the error in case of GLPK fit;
		 ErrorDriftRadiusbis[i]=ErrorDriftRadius[i]= 1. ;
		} else{
		// the following error is conventional for the chi**2 type of fit;
		 ErrorDriftRadiusbis[i]=ErrorDriftRadius[i]= 0.5 ;
		}
	}
	// the Mvd Strips hit
	for(j=0, i = fnMvdPixelHitsinTrack[ncand]; j< fnMvdStripHitsinTrack[ncand]; j++){
		k=fListMvdStripHitsinTrack[ncand][j];
		ZEDbis[i] = ZED[i] = fZMvdStrip[k];
		S[i] = atan2( fYMvdStrip[k]-fOy[ncand],fXMvdStrip[k]-fOx[ncand]);
		if(S[i]<0.) S[i] +=2.*PI;
		Sbis[i] = S[i] ;
		// DriftRadius is set conventionally at -1, for later use in the SZ fit;
		// the error on the point used in the fit is ErrorDriftRadius and this
		// is overestimated to be  1cm.
		DriftRadiusbis[i]=DriftRadius[i]=-1.;
		if(YesGLPKfitSZ){
		// the following is the error in case of GLPK fit;
		 ErrorDriftRadiusbis[i]=ErrorDriftRadius[i]= 1. ;
		} else {
		// the following error is conventional for the chi**2 type of fit;
		 ErrorDriftRadiusbis[i]=ErrorDriftRadius[i]= 0.5 ;
		}
		i ++;
	}

	// the SciTil hit ( when they are 2, the fS_SciTilHitsinTrack is already a mean
	// of the two; then consider only 1 SciTil hit, the first, and make an average
	// of the two Z positions).

	i = fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand];
	if(fnSciTilHitsinTrack[ncand]==2){
		ZED[i]=0.5*(fposizSciTil[fListSciTilHitsinTrack[ncand][0]][2]+
			fposizSciTil[fListSciTilHitsinTrack[ncand][1]][2]);
		S[i] = fS_SciTilHitsinTrack[ncand][0];
		// DriftRadius is set conventionally at -2, for later use in the SZ fit;
		// the error on the point used in the fit is ErrorDriftRadius and this
		// is overestimated to be DIMENSIONSCITIL/2.
		DriftRadius[i]=-2.;
		ErrorDriftRadius[i]= DIMENSIONSCITIL/2.; ;
	} else if (fnSciTilHitsinTrack[ncand]==1){
		ZED[i]=fposizSciTil[fListSciTilHitsinTrack[ncand][0]][2];
		S[i] = fS_SciTilHitsinTrack[ncand][0];
		// DriftRadius is set conventionally at -2, for later use in the SZ fit;
		// the error on the point used in the fit is ErrorDriftRadius and this
		// is overestimated to be DIMENSIONSCITIL/2.
		DriftRadius[i]=-2.;
		ErrorDriftRadius[i]= DIMENSIONSCITIL/2.; ;
	}



	// the Skew Stt hits
	// the info on TemporaryS, TemporaryZ, TemporaryZDrift, TemporaryZErrorafterTilt  as follows :
	// if   i  is the (original) Skew Hit number, and ii (0 or 1) is
	// the solution according to the calculateintersections method ,
	// then the infos are stored in the   i + ii*MAXSTTHITS location;

	for(j=0;j<fnSttSkewHitsinTrack[ncand]; j++){

		kall = fnMvdPixelHitsinTrack[ncand]+
			fnMvdStripHitsinTrack[ncand]+j;
		if( fnSciTilHitsinTrack[ncand] ==2 ){ // in this case only 1 SciTil hit
						// has been considered above;
		  i = fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+1+j;
		} else {  // this is the case of 1 or 0 SciTil hits in track; the
			//  (impossible?) case of > 2 SciTil hits has already been
			//  prevented early in PndTrkCTFindTrackInXY.
		  i = fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+
			fnSciTilHitsinTrack[ncand] +j ;
		}
		// here :		k = Skew Hit number (original notation),
		//  fListSttSkewHitsinTrackSolution[ncand][j] = number of the solution
		// (0 or 1) according to the  calculateintersections method;
		k=fListSttSkewHitsinTrack[ncand][j];
		location = k + fListSttSkewHitsinTrackSolution[ncand][j]*MAXSTTHITS;

		// calculate the quantities used for the SZ fit only.
		if(i<nhitsinfit){
			S[i] = TemporaryS[location];
			ZED[i]=TemporaryZ[location];
			DriftRadius[i]=TemporaryZDrift[location];
			ErrorDriftRadius[i]=2.*TemporaryZDrift[location];
		}
		//-----------
		ZEDbis[kall]=TemporaryZ[location];
		Sbis[kall]= TemporaryS[location];
		DriftRadiusbis[kall]=TemporaryZDrift[location];
		// overestimate the error on the Drift Radius used in the SZ  fit.
		if( fabs(TemporaryZDrift[location]) >1.e-10) {
		   ErrorDriftRadiusbis[kall]=TemporaryZErrorafterTilt[location];
		} else {
		   ErrorDriftRadiusbis[kall]=0.5;
		}


	}	//   end of  for(j=0;j<fnSttSkewHitsinTrack[ncand]; j++)


 return;
}
//----------end function PndTrkTracking::LoadSZetc_forSZfit



//--------------------------------  begin function   PndTrkTracking::MakeInclusionListStt


void PndTrkTracking::MakeInclusionListStt(
	Int_t nSttHit,
	Double_t info[][7]
	)
{

 int i,j;

// it needs to be initialized for each event !
 for(i=0;i<nSttHit;i++){
	fInclusionListStt[i] = true;
	SingleHitListStt[i] = true;
 }



//	fill the inclusion list for Stt, include only first hit for those straws with
//	multiple hits.


 for(i=0; i< nSttHit-1; i++){
	if( !fInclusionListStt[ i ] ) continue;
	for(j=i+1; j< nSttHit; j++){


		if(fInclusionListStt[ j ] &&
			fabs(info[i][0] - info[j][0])<1.e-20 &&
			fabs(info[i][1] - info[j][1])<1.e-20  )
		{
			SingleHitListStt[j]=fInclusionListStt[j]= false ;
		}
	} //  end of  for(j=i+1; j< Nhits;; j++)


 }   //   end of for(i=0; i< Nhits-1; i++)




}
//--------------------------------  end of function   PndTrkTracking::MakeInclusionListStt





//------------------------- begin of function  PndTrkTracking::MatchMvdHitsToSttTracks

void PndTrkTracking::MatchMvdHitsToSttTracks(
	Vec <bool>& keepit,
	Double_t delta,
	Double_t highqualitycut,
	Short_t nSttTrackCand,
	Double_t *FI0,
	Double_t *Fifirst,
	Vec <Short_t>& CHARGE,
	Short_t *nPixelHitsinTrack, // output
	Short_t ListPixelHitsinTrack[][MAXMVDPIXELHITSINTRACK], // output
	Short_t *nStripHitsinTrack, // output
	Short_t ListStripHitsinTrack[][MAXMVDSTRIPHITSINTRACK] // output
	)
{
	bool specialcase;

	Short_t i,
		jmvdhit;

	Double_t angle,
		anglemax,
		anglemin,
		dist;



 for(i=0; i<nSttTrackCand; i++){
	if( ! keepit[i] ) continue ;

	if( Fifirst[i] < -99998. ){  // case with Fifirst[i]=-99999.; in this
					// case there the circle is contained
					// in the Mvd region.
		anglemax = 2.*PI;
		anglemin = 0.;

	} else { // continuation of if( Fifirst[i] < -99998. )


	if(CHARGE[i]>0){	// track must rotate clockwise looking into the beam.
		anglemax = FI0[i];
		anglemin = Fifirst[i];
	} else {
		anglemin = FI0[i];
		anglemax = Fifirst[i];
	}
	if(anglemax < anglemin) anglemax += 2.*PI;
	if(anglemax < anglemin) anglemax=anglemin; // this is just to be super-sure.

	} // end of if( Fifirst[i] < -99998. )


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


 // first try attach the Pixels;

	nPixelHitsinTrack[i] = 0;
	for( jmvdhit=0; jmvdhit<fnMvdPixelHit; jmvdhit++){
		angle = atan2(fYMvdPixel[jmvdhit]-fOy[i],fXMvdPixel[jmvdhit]-fOx[i]);
		if(angle<0.) angle += 2.*PI;
		if( angle>anglemax){
			angle -= 2.*PI;
			if( angle>anglemax) angle = anglemax;

		} else if (angle<anglemin){
			angle += 2.*PI;
			if (angle<anglemin) angle = anglemin;
		}

		if(angle > anglemin && angle < anglemax)
		{
			dist=fabs( sqrt( (fOx[i]-fXMvdPixel[jmvdhit])* (fOx[i]-fXMvdPixel[jmvdhit])
			+(fOy[i]-fYMvdPixel[jmvdhit])*(fOy[i]-fYMvdPixel[jmvdhit]))-fR[i]);
			if(dist<delta)
			{
				ListPixelHitsinTrack[i][nPixelHitsinTrack[i]]=jmvdhit;
				nPixelHitsinTrack[i]++;
				if( nPixelHitsinTrack[i] == MAXMVDPIXELHITSINTRACK ) break;
			}
		}	// end of  if(angle > anglemin)
	 }  // end of  for( jmvdhit=0; jmvdhit<fnMvdPixelHits; jmvdhit++)


	// then try attach the Strips;
	nStripHitsinTrack[i] = 0;
	for( jmvdhit=0; jmvdhit<fnMvdStripHit; jmvdhit++){
		angle = atan2(fYMvdStrip[jmvdhit]-fOy[i],fXMvdStrip[jmvdhit]-fOx[i]);
		if(angle<0.) angle += 2.*PI;
		if( angle>anglemax){
			angle -= 2.*PI;
			if( angle>anglemax) angle = anglemax;
		} else if (angle<anglemin){
			angle += 2.*PI;
			if (angle<anglemin) angle = anglemin;
		}

		if(angle > anglemin && angle < anglemax)
		{
			dist=fabs( sqrt( (fOx[i]-fXMvdStrip[jmvdhit])* (fOx[i]-fXMvdStrip[jmvdhit])
			+(fOy[i]-fYMvdStrip[jmvdhit])*(fOy[i]-fYMvdStrip[jmvdhit]))-fR[i]);
			if(dist<delta)
			{
				ListStripHitsinTrack[i][nStripHitsinTrack[i]]=jmvdhit;
				nStripHitsinTrack[i]++;
				if( nStripHitsinTrack[i] == MAXMVDSTRIPHITSINTRACK ) break;
			}
		}	// end of  if(angle > anglemin)
	}  // end of  for( jmvdhit=0; jmvdhit<fnMvdStripHits; jmvdhit++)

 }	// end of for(i=0; i<nSttTrackCand; i++)


 return;
}

//------------------------- end of function  PndTrkTracking::MatchMvdHitsToSttTracks





//------------------------- begin of function  PndTrkTracking::MatchMvdHitsToSttTracksagain

void PndTrkTracking::MatchMvdHitsToSttTracksagain(
	Vec <bool>& keepit,
	Vec <bool>& Mvdhits,
	Double_t delta,
	Double_t highqualitycut,
	Short_t nSttTrackCand,
	Double_t *FI0,
	Double_t *Fifirst,
	Vec <Short_t>& CHARGE,

	Short_t *nPixelHitsinTrack, // output
	Short_t ListPixelHitsinTrack[][MAXMVDPIXELHITSINTRACK], // output
	Short_t *nStripHitsinTrack, // output
	Short_t ListStripHitsinTrack[][MAXMVDSTRIPHITSINTRACK] // output
	)
{
	bool
//		flaggo,
		specialcase,
	     downstream,
	     determined;

	Short_t j,
		itrack,
		ipix,
		istr,
		ndownstream,
		ntot,
		naddpix,
		naddstr,
		List[MAXMVDPIXELHITS+MAXMVDSTRIPHITS];

	Double_t angle,
		anglemax,
		anglemin,
		dist,
		oldtotal,
		oldtotal2,
		total,
		Dist,
		DIST[MAXMVDTRACKSPEREVENT+1];


  for(itrack=0; itrack<nSttTrackCand; itrack++){
	if( ! keepit[itrack] ) continue;
//	if( ! Mvdhits[itrack] ) continue;
	ntot=nPixelHitsinTrack[itrack]+nStripHitsinTrack[itrack];
	if( Fifirst[itrack] < -99998. ){  // case with Fifirst[i]=-99999.; in this
					// case the circle is contained
					// in the Mvd region.
		anglemax = 2.*PI;
		anglemin = 0.;
	} else {
		if(CHARGE[itrack]>0){
			// track must rotate clockwise looking into the beam.
			anglemax = FI0[itrack];
			anglemin = Fifirst[itrack];
		} else {
			anglemin = FI0[itrack];
			anglemax = Fifirst[itrack];
		}
		if(anglemax < anglemin) anglemax += 2.*PI;
		if(anglemax < anglemin) anglemax=anglemin;// just to be super-sure.

	} // end of if( Fifirst[itrack] < -99998. )


//  find if this track goes downstream or upstream

	ndownstream=0;
	for(j=0; j<nPixelHitsinTrack[itrack]; j++){
		if( fZMvdPixel[ListPixelHitsinTrack[itrack][j]]>0.){
			ndownstream++ ;
		}
	}
	for(j=0; j<nStripHitsinTrack[itrack]; j++){
		if( fZMvdPixel[ListStripHitsinTrack[itrack][j]]>0.){
			ndownstream++ ;
		}
	}
	if(ndownstream>ntot-ndownstream) downstream=true;
	else downstream=false;

//  loop over the Mvd Pixel and try to attach new Pixels to each candidate track

	naddpix=0;
	for(ipix=0; ipix<fnMvdPixelHit; ipix++){

/*
		flaggo=true;
		for(j=0; j<nPixelHitsinTrack[itrack]; j++){
		  if( ListPixelHitsinTrack[itrack][j]==ipix){
			flaggo=false;
			break;
		  }
		}
*/

//		if(flaggo){
		  angle = atan2(fYMvdPixel[ipix]-fOy[itrack],
			fXMvdPixel[ipix]-fOx[itrack]);
		  if(angle<0.) angle += 2.*PI;
		  if( angle>anglemax){
			angle -= 2.*PI;
			if( angle>anglemax) angle = anglemax;
		  } else if (angle<anglemin){
			angle += 2.*PI;
			if (angle<anglemin) angle = anglemin;
		  }
		  if(angle > anglemin && angle < anglemax){
			dist=fabs( sqrt(
			(fOx[itrack]-fXMvdPixel[ipix])*(fOx[itrack]-fXMvdPixel[ipix])
			+(fOy[itrack]-fYMvdPixel[ipix])*(fOy[itrack]-fYMvdPixel[ipix])
			) -fR[itrack]);
			if(dist<highqualitycut){
				List[naddpix]=ipix;
				naddpix++;
			}
		  }	// end of  if(angle > anglemin)
//		} // end of if(flaggo)

	}	//  end of   for(ipix=0; ipix<fnMvdPixelHit; ipix++)

	if(naddpix>0){
		if(naddpix>MAXMVDPIXELHITSINTRACK){
			// protection against strange tracks (and also from
			// out-of-bound indexing of arrays);
			naddpix=MAXMVDPIXELHITSINTRACK;
		}  // 
		for(j=0;j<naddpix;j++){
			ListPixelHitsinTrack[itrack][j]=List[j];
		}
		nPixelHitsinTrack[itrack] = naddpix;
	}	//  end of  if(naddpix>0)


	naddstr=0;
	for(istr=0; istr<fnMvdStripHit; istr++){
/*
		flaggo=true;
		for(j=0; j<nStripHitsinTrack[itrack]; j++){
		     if( ListStripHitsinTrack[itrack][j]==istr){
			flaggo=false;
			break;
		     }
		}
*/
//		if(flaggo){
		  angle = atan2(fYMvdStrip[istr]-fOy[itrack],
				fXMvdStrip[istr]-fOx[itrack]);
		  if(angle<0.) angle += 2.*PI;
		  if( angle>anglemax){
			angle -= 2.*PI;
			if( angle>anglemax) angle = anglemax;
		  } else if (angle<anglemin){
			angle += 2.*PI;
			if (angle<anglemin) angle = anglemin;
		  }
		  if(angle > anglemin && angle < anglemax){
			dist=fabs( sqrt(
			  (fOx[itrack]-fXMvdStrip[istr])*(fOx[itrack]-fXMvdStrip[istr])
			  +(fOy[itrack]-fYMvdStrip[istr])*(fOy[itrack]-fYMvdStrip[istr])
				) -fR[itrack]);

			if(dist<highqualitycut){
				List[naddstr]=istr;
				naddstr++;
			}
		  }	//  end of  if(angle > anglemin)

//		} // end of if(flaggo)

	}	//  end of   for(istr=0; istr<fnMvdStripHit; istr++)
	if(naddstr>0){
		if(naddstr>MAXMVDSTRIPHITSINTRACK){
			// protection against strange tracks (and also from
			// out-of-bound indexing of arrays);
			naddstr=MAXMVDSTRIPHITSINTRACK;
		}  // 
		for(j=0;j<naddstr;j++){
			ListStripHitsinTrack[itrack][j]=List[j];
		}
		nStripHitsinTrack[itrack] = naddstr;
	}	//  end of  if(naddpix>0)


	if(nPixelHitsinTrack[itrack]+nStripHitsinTrack[itrack]>0)
		Mvdhits[itrack]=true;

 }	// end of for(itrack=0; itrack<nSttTrackCand; itrack++)

 return;
}

//------------------- end function  PndTrkTracking::MatchMvdHitsToSttTracksagain




//------------------------- begin of function  PndTrkTracking::MatchMvdHitsToSttTracks2

void PndTrkTracking::MatchMvdHitsToSttTracks2(
	bool *keepit,
	Double_t delta,
	Double_t highqualitycut,
	Short_t nSttTrackCand,
	Double_t *FI0,
	Double_t *Fifirst,
	Short_t *CHARGE,
	Short_t *nPixelHitsinTrack, // output
	Short_t ListPixelHitsinTrack[][MAXMVDPIXELHITSINTRACK], // output
	Short_t *nStripHitsinTrack, // output
	Short_t ListStripHitsinTrack[][MAXMVDSTRIPHITSINTRACK] // output
	)
{
	bool specialcase;

	Short_t i,j,j1, j2, imvdcand, jmvdhit, ncont,
		chosenmix,
		chosenmix2,
		ngoodmix,
		oldN,
		nn[MAXMVDTRACKSPEREVENT+2],
		nHighQuality[MAXMVDTRACKSPEREVENT+2],
		List[MAXMVDTRACKSPEREVENT+2][MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK],
		ListType[MAXMVDTRACKSPEREVENT+2][MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK];

	Double_t angle,
		anglemax,
		anglemin,
		dist,
		oldtotal,
		oldtotal2,
		total,
		Dist,
		DIST[MAXMVDTRACKSPEREVENT+1];



 for(i=0; i<nSttTrackCand; i++){
	if( ! keepit[i] ) continue ;

	if( Fifirst[i] < -99998. ){  // case with Fifirst[i]=-99999.; in this
					// case there the circle is contained
					// in the Mvd region.
		anglemax = 2.*PI;
		anglemin = 0.;

	} else { // continuation of if( Fifirst[i] < -99998. )


	if(CHARGE[i]>0){	// track must rotate clockwise looking into the beam.
		anglemax = FI0[i];
		anglemin = Fifirst[i];
	} else {
		anglemin = FI0[i];
		anglemax = Fifirst[i];
	}
	if(anglemax < anglemin) anglemax += 2.*PI;
	if(anglemax < anglemin) anglemax=anglemin; // this is just to be super-sure.

	} // end of if( Fifirst[i] < -99998. )


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

	ngoodmix=0;
	nn[0]=0;
 for( imvdcand=0; imvdcand<fnMvdTrackCand; imvdcand++){
	Dist = 0.;
	ncont=0;
	nn[ngoodmix]=0;
	nHighQuality[ngoodmix]=0;
	for( jmvdhit=0; jmvdhit<fnHitMvdTrackCand[imvdcand]; jmvdhit++){

		if(fListHitTypeMvdTrackCand[imvdcand][jmvdhit]==
		    FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)){
			ncont++;
			angle = atan2(
			fYMvdPixel[fListHitMvdTrackCand[imvdcand][jmvdhit]]-fOy[i],
			fXMvdPixel[fListHitMvdTrackCand[imvdcand][jmvdhit]]-fOx[i]
							);
			if(angle<0.) angle += 2.*PI;

			if( angle>anglemax){
				angle -= 2.*PI;
				if( angle>anglemax) angle = anglemax;
			} else if (angle<anglemin){
				angle += 2.*PI;
				if (angle<anglemin) angle = anglemin;
			}
			if(angle > anglemin && angle < anglemax)
			{
				dist=fabs( sqrt(
				 (fOx[i]-fXMvdPixel[fListHitMvdTrackCand[imvdcand][jmvdhit]])*
				 (fOx[i]-fXMvdPixel[fListHitMvdTrackCand[imvdcand][jmvdhit]])
				+(fOy[i]-fYMvdPixel[fListHitMvdTrackCand[imvdcand][jmvdhit]])*
				 (fOy[i]-fYMvdPixel[fListHitMvdTrackCand[imvdcand][jmvdhit]]))-fR[i]);
				if(dist<delta)
				{
				     List[ngoodmix][nn[ngoodmix]]=
					fListHitMvdTrackCand[imvdcand][jmvdhit];
				     ListType[ngoodmix][nn[ngoodmix]]=
					FairRootManager::Instance()->GetBranchId(fMvdPixelBranch);
				     Dist += dist;
				     if( dist<highqualitycut) nHighQuality[ngoodmix]++;
				     nn[ngoodmix]++;
				}
			}	// end of  if(angle > anglemin)

		} else {// at this point this is a Strip hit; already made sure
			// earlier in the code that there is no third possibility.

			ncont++;
			angle = atan2(
			fYMvdStrip[fListHitMvdTrackCand[imvdcand][jmvdhit]]-fOy[i],
			fXMvdStrip[fListHitMvdTrackCand[imvdcand][jmvdhit]]-fOx[i]
							);
			if(angle<0.) angle += 2.*PI;

			if( angle>anglemax){
				angle -= 2.*PI;
				if( angle>anglemax) angle = anglemax;
			} else if (angle<anglemin){
				angle += 2.*PI;
				if (angle<anglemin) angle = anglemin;
			}
			if(angle > anglemin && angle < anglemax){
				dist=fabs( sqrt(
			 (fOx[i]-fXMvdStrip[fListHitMvdTrackCand[imvdcand][jmvdhit]])*
			 (fOx[i]-fXMvdStrip[fListHitMvdTrackCand[imvdcand][jmvdhit]])
			 +(fOy[i]-fYMvdStrip[fListHitMvdTrackCand[imvdcand][jmvdhit]])*
			 (fOy[i]-fYMvdStrip[fListHitMvdTrackCand[imvdcand][jmvdhit]])) -fR[i]);
				if(dist<delta)
				{
				   List[ngoodmix][nn[ngoodmix]]=
					fListHitMvdTrackCand[imvdcand][jmvdhit];
				   ListType[ngoodmix][nn[ngoodmix]]=
				    FairRootManager::Instance()->GetBranchId(fMvdStripBranch);
				   Dist += dist;
				   if( dist<highqualitycut) nHighQuality[ngoodmix]++;

				   nn[ngoodmix]++;
				}
			}	// end of   if(angle > anglemin)
		} // end of    if(fListHitTypeMvdTrackCand[imvdcand][jmvdhit]

	}	// end of   for( jmvdhit=0; jmvdhit<fnHitMvdTrackCand[imvdcand];


	if( nn[ngoodmix]>0) {
		DIST[ngoodmix]=Dist/nn[ngoodmix];
		ngoodmix++;


//--------- stampaggi
if(istampa>=3 ){cout<<"\tquesto Mvd candidato (n. ngoodmix = "<<ngoodmix-1<<
	") passa con i seguenti hits :"<<endl;

	for(int icc=0; icc<nn[ngoodmix-1]; icc++){
		if(ListType[ngoodmix-1][icc]==FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) {
			cout<<"\tPixel hit n. "<<List[ngoodmix-1][icc]<<endl;
		} else if(ListType[ngoodmix-1][icc]==
			FairRootManager::Instance()->GetBranchId(fMvdStripBranch)){
			cout<<"\tStrip hit n. "<<List[ngoodmix-1][icc]<<endl;
		} else{
			cout<<"\tNoise  (?) , hit type "<<ListType[ngoodmix-1][icc]<<endl;
		}
	}
	cout<<endl;
}
//------------fine stampaggi



	}
 }	// end of for( imvdcand=0;imvdcand<fnMvdTrackCand;imvdcand++)


//-------  now use the Mvd which are in no Mvd Track Candidate


//------- first, the DS (downstream) Mvd hits
		nn[ngoodmix]=0;
		DIST[ngoodmix] = 0.;
		nHighQuality[ngoodmix]=0;
		for( jmvdhit=0; jmvdhit<fnMvdDSPixelHitNotTrackCand; jmvdhit++){

			angle = atan2(
				fYMvdPixel[fListMvdDSPixelHitNotTrackCand[jmvdhit]]-fOy[i],
				fXMvdPixel[fListMvdDSPixelHitNotTrackCand[jmvdhit]]-fOx[i]
					);
			if(angle<0.) angle += 2.*PI;
			if( angle>anglemax){
				angle -= 2.*PI;
				if( angle>anglemax) angle = anglemax;
			} else if (angle<anglemin){
				angle += 2.*PI;
				if (angle<anglemin) angle = anglemin;
			}
			if(angle > anglemin && angle < anglemax){
				dist=fabs( sqrt(
				 (fOx[i]-fXMvdPixel[fListMvdDSPixelHitNotTrackCand[jmvdhit]])*
				 (fOx[i]-fXMvdPixel[fListMvdDSPixelHitNotTrackCand[jmvdhit]])
				+(fOy[i]-fYMvdPixel[fListMvdDSPixelHitNotTrackCand[jmvdhit]])*
				 (fOy[i]-fYMvdPixel[fListMvdDSPixelHitNotTrackCand[jmvdhit]])
						) -fR[i]);
				if(dist<delta)
				{
					List[ngoodmix][nn[ngoodmix]]=
					fListMvdDSPixelHitNotTrackCand[jmvdhit];
					ListType[ngoodmix][nn[ngoodmix]]=
					 FairRootManager::Instance()->GetBranchId(fMvdPixelBranch);
					DIST[ngoodmix] += dist;
					if( dist<highqualitycut) nHighQuality[ngoodmix]++;
						nn[ngoodmix]++;
				}
			}  //  end of     if(angle > anglemin )
		}	//  end  of for( jmvdhit=0; jmvdhit<fnMvdDSPixelHitNotTrackCand; jmvdhit++)

		for( jmvdhit=0; jmvdhit<fnMvdDSStripHitNotTrackCand; jmvdhit++){

			angle = atan2(
				fYMvdStrip[fListMvdDSStripHitNotTrackCand[jmvdhit]]-fOy[i],
				fXMvdStrip[fListMvdDSStripHitNotTrackCand[jmvdhit]]-fOx[i]
				      );
			if(angle<0.) angle += 2.*PI;
			if( angle>anglemax){
				angle -= 2.*PI;
				if( angle>anglemax) angle = anglemax;
			} else if (angle<anglemin){
				angle += 2.*PI;
				if (angle<anglemin) angle = anglemin;
			}
			if(angle > anglemin && angle < anglemax){

				dist=fabs( sqrt(
				 (fOx[i]-fXMvdStrip[fListMvdDSStripHitNotTrackCand[jmvdhit]])*
				 (fOx[i]-fXMvdStrip[fListMvdDSStripHitNotTrackCand[jmvdhit]])
				+(fOy[i]-fYMvdStrip[fListMvdDSStripHitNotTrackCand[jmvdhit]])*
				 (fOy[i]-fYMvdStrip[fListMvdDSStripHitNotTrackCand[jmvdhit]])
					) -fR[i]);


				if(dist<delta)
				{
					List[ngoodmix][nn[ngoodmix]]=
					fListMvdDSStripHitNotTrackCand[jmvdhit];
					ListType[ngoodmix][nn[ngoodmix]]=
					 FairRootManager::Instance()->GetBranchId(fMvdStripBranch);
					DIST[ngoodmix] += dist;
					if( dist<highqualitycut) nHighQuality[ngoodmix]++;
					nn[ngoodmix]++;
				}
			}  //  end of     if(angle > anglemin )

		}	//  end  of for( jmvdhit=0; jmvdhit<fnMvdDSStripHitNotTrackCand; jmvdhit++)

			if( nn[ngoodmix]>0) {
				DIST[ngoodmix] /= nn[ngoodmix];
				ngoodmix++;
//--------- stampaggi
if(istampa>=3){cout<<"\tevento n. "<<IVOLTE<<" questi Mvd ALONE DS hits passano  :\n"<<endl;

	for(int icc=0; icc<nn[ngoodmix-1]; icc++){
		if(ListType[ngoodmix-1][icc]==FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) {
			cout<<"\tDS Pixel hit n. "<<List[ngoodmix-1][icc]<<endl;
		} else if(ListType[ngoodmix-1][icc]==
			FairRootManager::Instance()->GetBranchId(fMvdStripBranch)){
			cout<<"\tDS Strip hit n. "<<List[ngoodmix-1][icc]<<endl;
		} else{
			cout<<"\tNoise (?) , hit tipo "<<ListType[ngoodmix-1][icc]<<endl;
		}
	}

}
//------------fine stampaggi
			}	// end of if( nn[ngoodmix]>0)






//--------  now the US (upstream) Mvd hits

		nn[ngoodmix]=0;
		DIST[ngoodmix] = 0.;
		nHighQuality[ngoodmix]=0;
		for( jmvdhit=0; jmvdhit<fnMvdUSPixelHitNotTrackCand; jmvdhit++){

			angle = atan2(
			fYMvdPixel[fListMvdUSPixelHitNotTrackCand[jmvdhit]]-fOy[i],
					fXMvdPixel[fListMvdUSPixelHitNotTrackCand[jmvdhit]]-fOx[i]
							);
			if(angle<0.) angle += 2.*PI;
			if( angle>anglemax){
				angle -= 2.*PI;
				if( angle>anglemax) angle = anglemax;
			} else if (angle<anglemin){
				angle += 2.*PI;
				if (angle<anglemin) angle = anglemin;
			}
			if(angle > anglemin && angle < anglemax){
				dist=fabs( sqrt(
				 (fOx[i]-fXMvdPixel[fListMvdUSPixelHitNotTrackCand[jmvdhit]])*
				 (fOx[i]-fXMvdPixel[fListMvdUSPixelHitNotTrackCand[jmvdhit]])
				+(fOy[i]-fYMvdPixel[fListMvdUSPixelHitNotTrackCand[jmvdhit]])*
				 (fOy[i]-fYMvdPixel[fListMvdUSPixelHitNotTrackCand[jmvdhit]])
						) -fR[i]);
				if(dist<delta)
				{
					List[ngoodmix][nn[ngoodmix]]=
						fListMvdUSPixelHitNotTrackCand[jmvdhit];
					ListType[ngoodmix][nn[ngoodmix]]=
					 FairRootManager::Instance()->GetBranchId(fMvdPixelBranch);
					DIST[ngoodmix] += dist;
					if( dist<highqualitycut) nHighQuality[ngoodmix]++;
						nn[ngoodmix]++;
				}
			}  //  end of     if(angle > anglemin )

		}	//  end  of for( jmvdhit=0; jmvdhit<fnMvdUSPixelHitNotTrackCand; jmvdhit++)

		for( jmvdhit=0; jmvdhit<fnMvdUSStripHitNotTrackCand; jmvdhit++){

			angle = atan2(
				fYMvdStrip[fListMvdUSStripHitNotTrackCand[jmvdhit]]-fOy[i],
				fXMvdStrip[fListMvdUSStripHitNotTrackCand[jmvdhit]]-fOx[i]
				      );
			if(angle<0.) angle += 2.*PI;
			if( angle>anglemax){
				angle -= 2.*PI;
				if( angle>anglemax) angle = anglemax;
			} else if (angle<anglemin){
				angle += 2.*PI;
				if (angle<anglemin) angle = anglemin;
			}
			if(angle > anglemin && angle < anglemax){

				dist=fabs( sqrt(
				 (fOx[i]-fXMvdStrip[fListMvdUSStripHitNotTrackCand[jmvdhit]])*
				 (fOx[i]-fXMvdStrip[fListMvdUSStripHitNotTrackCand[jmvdhit]])
				+(fOy[i]-fYMvdStrip[fListMvdUSStripHitNotTrackCand[jmvdhit]])*
				 (fOy[i]-fYMvdStrip[fListMvdUSStripHitNotTrackCand[jmvdhit]])
					) -fR[i]);


				if(dist<delta)
				{
					List[ngoodmix][nn[ngoodmix]]=
					fListMvdUSStripHitNotTrackCand[jmvdhit];
					ListType[ngoodmix][nn[ngoodmix]]=
					 FairRootManager::Instance()->GetBranchId(fMvdStripBranch);
					DIST[ngoodmix] += dist;
					if( dist<highqualitycut) nHighQuality[ngoodmix]++;
					nn[ngoodmix]++;
				}
			}  //  end of     if(angle > anglemin )

		}	//  end  of for( jmvdhit=0; jmvdhit<fnMvdUSStripHitNotTrackCand; jmvdhit++)

			if( nn[ngoodmix]>0) {
				DIST[ngoodmix] /= nn[ngoodmix];
				ngoodmix++;
//--------- stampaggi
if(istampa>=3){cout<<"\tevento n. "<<IVOLTE<<" questi Mvd ALONE US hits passano  :\n"<<endl;
	if( ngoodmix>0) {
	for(int icc=0; icc<nn[ngoodmix-1]; icc++){
		if(ListType[ngoodmix-1][icc]==FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)) {
			cout<<"\tUS Pixel hit n. "<<List[ngoodmix-1][icc]<<endl;
		} else if(ListType[ngoodmix-1][icc]==
			FairRootManager::Instance()->GetBranchId(fMvdStripBranch)){
			cout<<"\tUS Strip hit n. "<<List[ngoodmix-1][icc]<<endl;
		} else{
			cout<<"\tNoise (?) , hit tipo "<<ListType[ngoodmix-1][icc]<<endl;
		}
	}
	}

}
//------------fine stampaggi
			}	// end of if( nn[ngoodmix]>0)




//-------  end of using the Mvd which are in no Mvd Track Candidate


if(istampa>=3 ){cout<<"da PndTrkTracking : appena prima arbitration, IVOLTE = "<<
IVOLTE<<", Stt track cand = "<<i<<", ngoodmix = "<<ngoodmix<<endl;}

		if( ngoodmix==1){
			chosenmix=0;
			chosenmix2=0;
		} else if( ngoodmix>1) {
//--- here the arbitration if there are more than 1 Stt+Mvd hit combination for a given SttTrackCand
			oldtotal = DIST[0];
			oldtotal2 = DIST[0];
			oldN = nHighQuality[0];
//			oldtotal /= nTotali[0];
if(istampa>=3 ){cout<<"da PndTrkTracking : goodmix n. 0, total distance (che e' = total distance2) = "<<oldtotal
				<<", e nHighQuality = "<<nHighQuality[0]<<endl;}
			chosenmix=0;
			chosenmix2=0;
			for(j1=1; j1<ngoodmix;j1++){
				total = DIST[j1];
if(istampa>=3){cout<<"da PndTrkTracking :\t goodmix n. "<<j1<<", total distance "<<total
					<<", e nHighQuality = "<<nHighQuality[j1]<<endl;}
				if(oldN<nHighQuality[j1]){
					oldN=nHighQuality[j1];
					chosenmix2=j1;
				} else if (oldN==nHighQuality[j1]){
					if(total<oldtotal2){
						chosenmix2=j1;
						oldtotal2=total;
					}
				}
				if(total<oldtotal){
					oldtotal=total;
					chosenmix=j1;
				}
			}
		}	// end of  if( ngoodmix==1)
//--- end of arbitration
if(istampa>=3 ){cout<<"da PndTrkTracking : fine arbitration, IVOLTE = "<<
IVOLTE<<", Stt track cand = "<<i<<endl;}





	nPixelHitsinTrack[i]=0;
	nStripHitsinTrack[i]=0;
	if( ngoodmix>0){
		chosenmix=chosenmix2;
		for(j=0;j<nn[chosenmix];j++){
			if( ListType[chosenmix][j]==
				FairRootManager::Instance()->GetBranchId(fMvdPixelBranch) ){
				ListPixelHitsinTrack[i]
				   [nPixelHitsinTrack[i]]=List[chosenmix][j];
				nPixelHitsinTrack[i]++;
			} else if(
			     ListType[chosenmix][j]==
			     FairRootManager::Instance()->GetBranchId(fMvdStripBranch)){
				ListStripHitsinTrack[i]
				   [nStripHitsinTrack[i]]=List[chosenmix][j];
				nStripHitsinTrack[i]++;
			}
		}
	}	// end of if( ngoodmix>0)


 }	// end of for(i=0; i<nSttTrackCand; i++)


 return;
}

//------------------------- end of function  PndTrkTracking::MatchMvdHitsToSttTracks2


//------begin function PndTrkTracking::OrderingConformal_Loading_ListTrackCandHit
void PndTrkTracking::OrderingConformal_Loading_ListTrackCandHit(
	Vec <bool>& keepit,
	Short_t ncand,
	Double_t info[][7],
	Double_t Trajectory_Start[][2],
	Vec <Short_t>& CHARGE,
	Double_t SchosenSkew[][MAXSTTHITS]
				)
{

	Short_t	i,
			j,
			ipar,
			iskew;


//     ordering all the hits belonging to the candidate track, by increasing fR;
//     forming the new track with Mvd+Stt hits


	// arrays used to store temporarily the info of Mvd hits to be ordered.
	Int_t ListHits[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK],
		ListHits2[MAXSTTHITSINTRACK];
	Double_t XY[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK][2],
		XY2[MAXSTTHITSINTRACK][2];

//	for(ncand=FirstCandidate; ncand< LastCandidate; ncand++){

		if(!keepit[ncand]) return;
		fnTrackCandHit[ncand] =fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand]+
					fnMvdPixelHitsinTrack[ncand]+
					fnMvdStripHitsinTrack[ncand];
		// adding the Mvd hits (Pixel and Strips)

		if( fnMvdPixelHitsinTrack[ncand]+
		    fnMvdStripHitsinTrack[ncand] >0){
			for(i=0; i< fnMvdPixelHitsinTrack[ncand]; i++){
			  XY[i][0] = fXMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ];
			  XY[i][1] = fYMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ];
			  ListHits[i] = fListMvdPixelHitsinTrack[ncand][i];
			}
			for(i=0; i< fnMvdStripHitsinTrack[ncand]; i++){
			  XY[i+fnMvdPixelHitsinTrack[ncand]][0] =
			    fXMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ];
			  XY[i+fnMvdPixelHitsinTrack[ncand]][1] =
			    fYMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ];
			  // to distinguish between Pixels and Strips, add a number
			  // to the original Strip hit number.
			  ListHits[i+fnMvdPixelHitsinTrack[ncand]] =
				fListMvdStripHitsinTrack[ncand][i]+
				(MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*10 ;
			}

			//  ordering the Mvd Hits
			OrderingUsingConformal(
			   fOx[ncand],
			   fOy[ncand],
			   &Trajectory_Start[ncand][0],
			   fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand],
			   XY, // XY[*][0] = X position, XY[*][0] = Y position.
			   CHARGE[ncand],  // input
			   ListHits  // output
						);
			//  constructing the ordered new Track  Candidate now
			for(i=0; i< fnMvdPixelHitsinTrack[ncand]+
					fnMvdStripHitsinTrack[ncand]; i++){
				if(ListHits[i]<(MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*10){//Pixel.
					fListTrackCandHit[ncand][i] = ListHits[i];
					fListTrackCandHitType[ncand][i] = 0;
				} else { // Strip hits.
					fListTrackCandHit[ncand][i] = ListHits[i]-
						(MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*10;
					fListTrackCandHitType[ncand][i] = 1;
				}
			}  // end of  for(i=0; i< fnMvdPixelHitsinTrack[ncand]+

		}	// end of  if( fnMvdPixelHitsinTrack[ncand]+




		// construction of the second part of the ordered new Track  Candidate

		if( fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand] >0){

			for(i=0; i<fnSttParHitsinTrack[ncand]; i++){
			 XY2[i][0] = info[ fListSttParHitsinTrack[ncand][i] ][0];
			 XY2[i][1] = info[ fListSttParHitsinTrack[ncand][i] ][1];
			 ListHits2[i] = fListSttParHitsinTrack[ncand][i];


			}

			for(i=0; i<fnSttSkewHitsinTrack[ncand]; i++){
			 j = i+fnSttParHitsinTrack[ncand];
			 XY2[j][0] = fOx[ncand]+fR[ncand]*
				cos(SchosenSkew[ncand][fListSttSkewHitsinTrack[ncand][i]]);
			 XY2[j][1] = fOy[ncand]+fR[ncand]*
				sin(SchosenSkew[ncand][fListSttSkewHitsinTrack[ncand][i]]);
			 ListHits2[j] = fListSttSkewHitsinTrack[ncand][i]+
					MAXSTTHITS*10; // in order to distinguish
							//  the Skew hits.
			}


		//  ordering the Stt Hits

			OrderingUsingConformal(
			 fOx[ncand],
			 fOy[ncand],
			 &Trajectory_Start[ncand][0],
			 fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand],
			 XY2, // XY2[*][0] = X position, XY2[*][0] = Y position.
			 CHARGE[ncand],  // input
			 ListHits2  // output
					);

			for(j=0,ipar=0,iskew=0;
				j< fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand];j++){
			  i = j+fnMvdPixelHitsinTrack[ncand]+ fnMvdStripHitsinTrack[ncand];
			  if(ListHits2[j]<MAXSTTHITS*10){  // parallel Stt hit.
			    fListTrackCandHit[ncand][i] = ListHits2[j];
			    fListTrackCandHitType[ncand][i] = 2;
			    fListSttParHitsinTrack[ncand][ipar]=ListHits2[j];
			    ipar++;
			  } else {  // skew Stt hit.
			    fListTrackCandHit[ncand][i] = ListHits2[j]-MAXSTTHITS*10;
			    fListTrackCandHitType[ncand][i] = 3;
			    fListSttSkewHitsinTrack[ncand][iskew]=ListHits2[j]-MAXSTTHITS*10;
			    iskew++;
			  }  // end of  if(ListHits2[j]<MAXSTTHITS*10)
			} // end of for(j=0,ipar=0,iskew=0; ....

		}	// end of  if( fnSttParHitsinTrack[ncand]+


if(istampa>=3) for(int ica=0; ica<fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand]+
	fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand]; ica++){
	cout<<"from PndTrkTracking, hit n. "<<fListTrackCandHit[ncand][ica]<<", hit type "
	<<fListTrackCandHitType[ncand][ica]<<endl;
}


//	} //   end of  for(ncand=FirstCandidate; ncand< LastCandidate; ncand++)


	return;
}

//--------end function PndTrkTracking::OrderingConformal_Loading_ListTrackCandHit

//----------begin of function PndTrkTracking::Ordering_Loading_ListTrackCandHit

void PndTrkTracking::Ordering_Loading_ListTrackCandHit(
	Vec<bool>& keepit,
	Short_t FirstCandidate,
	Short_t LastCandidate,
	Double_t info[][7],
	Double_t Trajectory_Start[][2],
	Vec <Short_t>& CHARGE,
	Double_t SchosenSkew[][MAXSTTHITS]
				)
{
	Short_t ncand;
	for(ncand=FirstCandidate; ncand< LastCandidate; ncand++){

		// for small radius trajectory better the ordering with conformal.
		if( fR[ncand]< RSTRAWDETECTORMAX/2.){
			OrderingConformal_Loading_ListTrackCandHit(
				keepit,
				ncand,
				info,
				Trajectory_Start,
				CHARGE,
				SchosenSkew
						);
		} else { // otherwise it is better distance from (0,0) method.
			OrderingR_Loading_ListTrackCandHit(
				keepit,
				ncand,
				info
						);
		}

	} //   end of  for(ncand=FirstCandidate; ncand< LastCandidate; ncand++)


	return;
}
//----------end of function PndTrkTracking::Ordering_Loading_ListTrackCandHit



//----------begin of function PndTrkTracking::OrderingR_Loading_ListTrackCandHit

void PndTrkTracking::OrderingR_Loading_ListTrackCandHit(
	Vec <bool>& keepit,
	Short_t ncand,
	Double_t info[][7]
	)
{

 Short_t	dime,
		i,
		j,
		ipar,
		iskew;

 PndTrkMergeSort MergeSort;

//     ordering all the hits belonging to the candidate track, by increasing fR;
//     forming the new track with Mvd+Stt hits


		if(!keepit[ncand]) return;
		fnTrackCandHit[ncand] =fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand]+
					fnMvdPixelHitsinTrack[ncand]+
					fnMvdStripHitsinTrack[ncand];
		dime = fnMvdPixelHitsinTrack[ncand]+fnMvdStripHitsinTrack[ncand];
		if(dime == 0 ) dime =1;	// just for protection;
		Short_t tempmvdindex[dime ],
			 tempmvdtype[dime ];
		Int_t auxIndex[dime ];
		Double_t auxR[dime ];
		// adding the Mvd hits (Pixel and Strips)
		for(i=0; i< fnMvdPixelHitsinTrack[ncand]; i++){
			auxR[i] =
			 fXMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ]*
			 fXMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ]+
			 fYMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ]*
			 fYMvdPixel[ fListMvdPixelHitsinTrack[ncand][i] ];
			tempmvdindex[i]=fListMvdPixelHitsinTrack[ncand][i];
			tempmvdtype[i]=0;
			auxIndex[i] = i;
		}
		for(i=0; i< fnMvdStripHitsinTrack[ncand]; i++){
			auxR[i+fnMvdPixelHitsinTrack[ncand]] =
			 fXMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ]*
			 fXMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ]+
			 fYMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ]*
			 fYMvdStrip[ fListMvdStripHitsinTrack[ncand][i] ];
			tempmvdindex[i+fnMvdPixelHitsinTrack[ncand]]=
				fListMvdStripHitsinTrack[ncand][i];
			tempmvdtype[i+fnMvdPixelHitsinTrack[ncand]]=1;
			auxIndex[i+fnMvdPixelHitsinTrack[ncand]]=
				i+fnMvdPixelHitsinTrack[ncand];
		}

		//  ordering the Mvd Hits
		if( fnMvdPixelHitsinTrack[ncand]+
		    fnMvdStripHitsinTrack[ncand] >0){
			MergeSort.Merge_Sort( fnMvdPixelHitsinTrack[ncand]+
		                    fnMvdStripHitsinTrack[ncand],
				    auxR, auxIndex);

		//  constructing the first part of the ordered new Track  Candidate
			for(i=0; i< fnMvdPixelHitsinTrack[ncand]+
					fnMvdStripHitsinTrack[ncand]; i++){
				fListTrackCandHit[ncand][i] = tempmvdindex[ auxIndex[i] ];
				fListTrackCandHitType[ncand][i] = tempmvdtype[ auxIndex[i] ];

			}
		}	// end of  if( fnMvdPixelHitsinTrack[ncand]+

		// construction of the second part of the ordered new Track  Candidate

		dime = fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand];
		if(dime == 0) dime =1;	// only for protection;
		Short_t tempmvdindex2[dime ],
			 tempmvdtype2[dime ];
		Int_t auxIndex2[dime ];
		Double_t auxR2[dime ];

		for(i=0; i<fnSttParHitsinTrack[ncand]; i++){
			auxR2[i] =
			 info[ fListSttParHitsinTrack[ncand][i] ][0]*
			 info[ fListSttParHitsinTrack[ncand][i] ][0]+
			 info[ fListSttParHitsinTrack[ncand][i] ][1]*
			 info[ fListSttParHitsinTrack[ncand][i] ][1];
			tempmvdindex2[i]=fListSttParHitsinTrack[ncand][i];
			tempmvdtype2[i]=2;
			auxIndex2[i] = i;
		}
		for(i=0; i<fnSttSkewHitsinTrack[ncand]; i++){
			j = i+fnSttParHitsinTrack[ncand];
			auxR2[j] =
			 info[ fListSttSkewHitsinTrack[ncand][i] ][0]*
			 info[ fListSttSkewHitsinTrack[ncand][i] ][0]+
			 info[ fListSttSkewHitsinTrack[ncand][i] ][1]*
			 info[ fListSttSkewHitsinTrack[ncand][i] ][1];
			tempmvdindex2[j]=fListSttSkewHitsinTrack[ncand][i];
			tempmvdtype2[j]=3;
			auxIndex2[j] = j;
		}


		//  ordering the Stt Hits
		if( fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand] >0){
			MergeSort.Merge_Sort( fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand],
				    auxR2, auxIndex2);

			for(j=0,ipar=0,iskew=0;
				  j< fnSttParHitsinTrack[ncand]+fnSttSkewHitsinTrack[ncand];j++){
				i = j+fnMvdPixelHitsinTrack[ncand]+
				    fnMvdStripHitsinTrack[ncand];
				fListTrackCandHit[ncand][i] = tempmvdindex2[ auxIndex2[j] ];
				fListTrackCandHitType[ncand][i] = tempmvdtype2[ auxIndex2[j] ];

				if( fListTrackCandHitType[ncand][i]==2) {
				  fListSttParHitsinTrack[ncand][ipar]=tempmvdindex2[auxIndex2[j]];
				  ipar++;
				} else {
				  fListSttSkewHitsinTrack[ncand][iskew]=tempmvdindex2[auxIndex2[j]];
				  iskew++;
				}
			}
		}	// end of  if( fnSttParHitsinTrack[ncand]+






	return;
}

//----------end of function PndTrkTracking::OrderingR_Loading_ListTrackCandHit



//----------begin of function PndTrkTracking::OrderingSttSkewandSttParallel

 void PndTrkTracking::OrderingSttSkewandSttParallel(
			Double_t oX,
			Double_t oY,
			Double_t Rr,
			Short_t nSkewhit,
			Short_t *ListSkewHits,
			Double_t *SList, // this is rekated to the skew hits. IMPORTANT :
					// the index must be the ORIGINAL skew hit number,
					// therefore SList[Infoskew[ListSkewHits[*]]].
			Short_t  Charge,
			Short_t nParHits,
			Short_t *ListParHits,
			Double_t *U,
			Double_t *V,
			Short_t *BigList // this is the final ordered Parallel+Skew list;
				// already in NATIVE hit number.
				)
{

 Short_t
	i,
	j,
	tmp[nSkewhit+nParHits],
	tmpList[nSkewhit];

 Int_t	index[nSkewhit+nParHits];

 Double_t	aaa,
		b1,
		sign,
		aux[nSkewhit+nParHits];


 PndTrkMergeSort MergeSort;


//  here there is the ordering of the hits, under the assumption that the circumference
//  in XY goes through (0,0).
//  Moreover, the code before is supposed to have selected trajectories in XY with (fOx,fOy)
//  farther from (0,0) by > 0.9 * RminStrawDetector/2 and consequently fOx and fOy are not both 0.
//  The scheme for the ordering of the hit is as follows :
//  1)  order hits by increasing U of the conformal mapping; see Gianluigi's Logbook page 283;
//  2)  find the charge of the track by checking if it is closest to the center in XY
//	the first or the last of the ordered hits.



//   ordering of the hits

	aaa = atan2( oY, oX);  // atan2 defined between -PI and PI.

	// the following statement is necessary since for unknown reason the root interpreter
	// gives a weird error when using PI directly in the if statement below!!!!!!! I lost
	// 2 hours trying to figure this out!
	b1 = PI/4.;


	if(aaa>b1&&aaa<3.*b1|| (aaa>-3.*b1&&aaa<-b1)){  //  case #1 or #3;see Gianluigi's Logbook page 285.
		if( (aaa>b1&&aaa<3.*b1 && Charge == -1)||( aaa>-3.*b1&&aaa<-b1 && Charge == 1) )
				{  // for speeding up the ordering taking advantage
				    // that the parallel hits were earlier ordered and
				    //  apply the trick of multiplying by   -1.
			sign=-1.;
		} else {  //  normal calculation
			sign=1.;
		}

		for (j = 0 ; j< nParHits; j++){
			aux[j] = sign*U[j];
		}
		for (j = 0; j< nSkewhit; j++){
			// this is U in conformal space
			aux[j+nParHits]=sign*(oX + Rr*cos(SList[ListSkewHits[j]]))/
					(oX*oX+oY*oY+Rr*Rr + 2.*Rr*
					(oX*cos(SList[ListSkewHits[j]])
					+oY*sin(SList[ListSkewHits[j]])));
		}


	} else {    // use V as ordering variable
		    // [case 2 and 4 Gianluigi's Logbook page 285].

		if( ((aaa<=-3.*b1|| aaa>=3.*b1) && Charge == -1)
				|| ( -b1 <= aaa && aaa <= b1 && Charge == 1) ){
			sign=-1.;
		} else {
			sign=1.;
		}
		for (j = 0 ; j< nParHits; j++){
			aux[j] = sign*V[j];
		}
		for (j = 0; j< nSkewhit; j++){
			// this is V in conformal space.
			aux[j+nParHits]=sign*(oY + Rr*sin(SList[ListSkewHits[j]]))/
			(oX*oX+oY*oY+Rr*Rr + 2.*Rr*
				(oX*cos(SList[ListSkewHits[j]])
				+oY*sin(SList[ListSkewHits[j]])));
		}


	}  //  end of  if((aaa>b1&& ....


	for (j = 0 ; j< nParHits; j++){
		BigList[j]=ListParHits[j];
		index[j] = j;
	}
	for (j = 0; j< nSkewhit; j++){
		BigList[j+nParHits]=ListSkewHits[j];
		index[j+nParHits] = j+nParHits;
	}

	MergeSort.Merge_Sort( nSkewhit+nParHits, aux, index);


	for(i=0, j=0;i<nSkewhit+nParHits;i++){
		tmp[i]=BigList[index[i]];
		//  reorder the ListSkewHits also.
		if( index[i] >= nParHits ){
			tmpList[j] = ListSkewHits[index[i]-nParHits];
			j++;
		}
	}
	for(i=0;i<nSkewhit+nParHits;i++){
		BigList[i]= tmp[i];
	}
		//  reorder the ListSkewHits also.
	for(i=0;i<nSkewhit;i++){
		ListSkewHits[i]= tmpList[i];
	}


 return;

}
//----------end of function PndTrkTracking::OrderingSttSkewandSttParallel


//----------begin of function PndTrkTracking::OrderingUsingConformal

void   PndTrkTracking::OrderingUsingConformal(
	Double_t oX,
	Double_t oY,
	Double_t Traj_Sta[2],
	Int_t nHits,
	Double_t XY[][2],
	Short_t  Charge,  // input
	Int_t *ListHits
	)
{




      Short_t	dime,i,j;
      Double_t	aaa,
		bbb,
		ccc,
		b1,
		firstR2,
		lastR2;

      dime = nHits;
      if(nHits==0) dime =1;

      Short_t	tmp[dime];
      Double_t	aux[dime],
		U[dime],
		V[dime];

 PndTrkMergeSort MergeSort;



//  here there is the ordering of the hits, NOT under the assumption that the circumference
//  in XY goes through  Trajectory_Start.
//  Moreover, the code before is supposed to have selected trajectories in XY with (fOx,fOy)
//  farther from (0,0) by > 0.9 * RminStrawDetector/2 and consequently fOx and fOy are not both 0.
//  The scheme for the ordering of the hit is as follows :
//  1)  order hits by increasing U or V of the conformal mapping; see Gianluigi's Logbook page 283;
//  2)  find the charge of the track by checking if it is closest to the center in XY
//	the first or the last of the ordered hits.
//  3)  in case, invert the ordering of U, V and ListHits such that the first hits in the
//	list are always those closer to the Trajectory_Start.


//   ordering of the hits

	aaa = atan2( oY-Traj_Sta[1], oX-Traj_Sta[0]);  // atan2 defined between -PI and PI.

	// the following statement is necessary since for unknown reason the root interpreter
	// gives a weird error when using PI directly in the if statement below!!!!!!! I lost
	// 2 hours trying to figure this out!
	b1 = PI/4.;

	if((aaa>b1&&aaa<3.*b1) || (aaa>-3.*b1&&aaa<-b1)){//use U as ordering variable;
							//[case 1 or 3 Gianluigi's Logbook page 285].
		for (j = 0; j< nHits; j++){
			bbb = XY[j][0]-Traj_Sta[0];
			ccc = XY[j][1]-Traj_Sta[1];
			U[j]= bbb/(bbb*bbb+ccc*ccc);
		}
		MergeSort.Merge_Sort( nHits, U, ListHits);

		if((aaa>b1&&aaa<3.*b1)){  //  case #1;
			if( Charge == -1){
				// inverting the order of the hits.
				for(i=0;i<nHits;i++){
					tmp[i]=ListHits[nHits-1-i];
				}
				for(i=0;i<nHits;i++){
					ListHits[i]=tmp[i];
				}
			}
		} else{  //  case # 3.
			if(Charge == 1){
				// inverting the order of the hits.
				for(i=0;i<nHits;i++){
					tmp[i]=ListHits[nHits-1-i];
				}
				for(i=0;i<nHits;i++){
					ListHits[i]=tmp[i];
				}
			}// end of  if( Charge ==1)
		}// end of  if((aaa>b1&&aaa<3.*b1))

	} else { // use V as ordering variable [case 2 or 4 Gianluigi's Logbook page 285].
		for (j = 0; j< nHits; j++){
			bbb = XY[j][0]-Traj_Sta[0];
			ccc = XY[j][1]-Traj_Sta[1];
			V[j]= ccc/(bbb*bbb+ccc*ccc);
		}
		MergeSort.Merge_Sort( nHits, V, ListHits);

		if((aaa<=-3.*b1 || aaa>=3.*b1)){  //  case #2;
			if( Charge == -1){
				// inverting the order of the hits.
				for(i=0;i<nHits;i++){
					tmp[i]=ListHits[nHits-1-i];
				}
				for(i=0;i<nHits;i++){
					ListHits[i]=tmp[i];
				}
			}
		} else{  //  case # 4.
			if( Charge == 1){
				// inverting the order of the hits.
				for(i=0;i<nHits;i++){
					tmp[i]=ListHits[nHits-1-i];
				}
				for(i=0;i<nHits;i++){
					ListHits[i]=tmp[i];
				}
			}
		}

	} //  end of   if((aaa>b1&& ....



 return; 


}
//----------end of function PndTrkTracking::OrderingUsingConformal


//------------------ begin function  PndTrkTracking::RefitMvdStt

void  PndTrkTracking::RefitMvdStt(
	Short_t nCandHit,
	Short_t *ListCandHit,
	Short_t *ListCandHitType,
	Double_t info[][7],
	Double_t rotationangle, //  this is between 0. and 2*PI
	Double_t tv[2],
	Short_t iexcl,
	Double_t *pAlfa, // output of the fit
	Double_t *pBeta, // output of the fit
	Double_t *pGamma,// set at zero always for now
	bool *status    // fit status; true = successful
	)
{
	bool Type;

	Short_t	i,
		iparallel;


	Short_t	exitstatus;

	Double_t dist2,
		 mindis,
		 emme,
		 factor,
		 gamma,
		 qu,
		 ErrorStraw = 0.03,
		 ErrorMvd = 0.01,
Xconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK],
Yconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK],
DriftRadiusconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK],
ErrorDriftRadiusconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK];


// PndTrkGlpkFits fit;
// PndTrkLegendreFits fit;
 PndTrkChi2Fits fit;


	*status= false;
	factor=3.;
	mindis=0.5;

	for(i=0, iparallel=0; i<nCandHit && iparallel< MAXHITSINFIT; i++){



		if(i==iexcl) continue;
		if( ListCandHitType[i] == 0 ){	// mvd pixels
			//----- translate the little circumference in XY representing
			//  approximately the sensitive area into the conformal space
			//  circumference

			dist2 = (fXMvdPixel[ListCandHit[i]]-tv[0])*
				(fXMvdPixel[ListCandHit[i]]-tv[0])+
				(fYMvdPixel[ListCandHit[i]]-tv[1])*
				(fYMvdPixel[ListCandHit[i]]-tv[1]);
			if(dist2<mindis) continue;	// this is to exclude Mvd hits too close to the traslated
						// center of reference frame (given by tv[0], tv[1]). This
						// would cause the subsequent fit to fail.

			gamma = dist2 - ErrorMvd*ErrorMvd; // for Pixels
				// I assume the 'drift radius' to be  the max dimension of
				//  the Pixel
			Xconformal[iparallel] = (fXMvdPixel[ListCandHit[i]]-tv[0])/gamma;
			Yconformal[iparallel] = (fYMvdPixel[ListCandHit[i]]-tv[1])/gamma;
			DriftRadiusconformal[iparallel]=-1.;// only to signal later this is a Mvd hit.
			ErrorDriftRadiusconformal[iparallel]=factor*ErrorMvd/fabs(gamma);
			iparallel++;
		} else if ( ListCandHitType[i] == 1 ){	// mvd strips
			//----- translate the little circumference in XY representing
			//  approximately the sensitive area into the conformal space
			//  circumference
			dist2 = (fXMvdStrip[ListCandHit[i]]-tv[0])*
				(fXMvdStrip[ListCandHit[i]]-tv[0])+
				(fYMvdStrip[ListCandHit[i]]-tv[1])*
				(fYMvdStrip[ListCandHit[i]]-tv[1]);
			if(dist2<mindis) continue;	// this is to exclude Mvd hits too close to the traslated
						// center of reference frame (given by tv[0], tv[1]). This
						// would cause the subsequent fit to fail.

			gamma = dist2 - ErrorMvd*ErrorMvd;// for Strips also
				// I assume the 'drift radius' to be  the largest error of
				//  the Strips (which is in the X dimension)
			Xconformal[iparallel] = (fXMvdStrip[ListCandHit[i]]-tv[0])/gamma;
			Yconformal[iparallel] = (fYMvdStrip[ListCandHit[i]]-tv[1])/gamma;
			DriftRadiusconformal[iparallel]=-1.;// only to signal later this is a Mvd hit.
			ErrorDriftRadiusconformal[iparallel]=factor* ErrorMvd/fabs(gamma);
			iparallel++;
		} else if ( ListCandHitType[i] == 2 ){	// Stt parallel hit.


			dist2 = (info[ListCandHit[i]][0]-tv[0])*
				(info[ListCandHit[i]][0]-tv[0])+
				(info[ListCandHit[i]][1]-tv[1])*
				(info[ListCandHit[i]][1]-tv[1]);
			// this is to exclude hits too close to the traslated
			// center of reference frame (given by tv[0], tv[1]). This
			// would cause the subsequent fit to fail.
			if(dist2<mindis) continue;
			gamma = dist2 -
				info[ListCandHit[i]][3]*
				info[ListCandHit[i]][3];
			Xconformal[iparallel] = (info[ListCandHit[i]][0]-tv[0])/gamma;
			Yconformal[iparallel] = (info[ListCandHit[i]][1]-tv[1])/gamma;
			DriftRadiusconformal[iparallel]=info[ListCandHit[i]][3]/fabs(gamma);
			ErrorDriftRadiusconformal[iparallel]=factor*ErrorStraw/fabs(gamma);
			iparallel++;

		}
	}	// end of for(i=0, iparallel=0;

	if( nCandHit < 2)  return;

	exitstatus = fit.FitHelixCylinder(
					iparallel,
					Xconformal,
					Yconformal,
					DriftRadiusconformal,
					ErrorDriftRadiusconformal,
					rotationangle,  //  rotationangle, da mettere
					tv,	//  vertex in (X,Y) of this trajectory
					MAXHITSINFIT,  //  maximum n. of hits allowed in fast fit
					&emme,
					&qu,
					pAlfa,
					pBeta,
					pGamma,
					&Type,
					0, // istampa
					0 //  IVOLTE
					);


	//  existatus > 0, Type= true --> fit ok, it is a Circle in XY; 
	//  existatus > 0, Type= false --> fit ok, it is a Straigh Line in XY; 
	//  existatus < 0, fit failed;  Type was set to false in this case;
	if( exitstatus > 0 && Type)	*status=true;
	return;
}


//------------------ end function  PndTrkTracking::RefitMvdStt

ClassImp(PndTrkTracking)