[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerFitter.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  Track fitter                                                             //
//                                                                           //
//  Authors:                                                                 //
//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
//    Markus Fasel <M.Fasel@gsi.de>                                          //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TROOT.h>
#include <TLinearFitter.h>
#include <TMath.h>
#include <TTreeStream.h>

#include "AliLog.h"

#include "AliTRDseed.h"
#include "AliTRDcalibDB.h"
#include "AliTRDcluster.h"
#include "AliTRDReconstructor.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackerFitter.h"

#define DEBUG

ClassImp(AliTRDtrackerFitter)

//_________________________________________________________________________
AliTRDtrackerFitter::AliTRDtrackerFitter()
  :TObject()
  ,fFitterTC(0x0)
  ,fFitterT2(0x0)
  ,fRieman1(0x0)
  ,fRieman2(0x0)
  ,fChi2TR(0)
  ,fChi2TC(0)
  ,fCR(0)
  ,fCC(0)
  ,fDca(0)
  ,fDzmf(0)
  ,fZmf(0)
  ,fNlayers(0)
  ,fDebugStream(0x0)
{
  //
  // Constructor
  //

	fRieman1   = new AliRieman(1000);
	fRieman2  = new AliRieman(1000);
	fFitterTC = new TLinearFitter(2,"hyp2");
	fFitterT2 = new TLinearFitter(4,"hyp4");
	fFitterTC->StoreData(kTRUE);
	fFitterT2->StoreData(kTRUE);
}

//_________________________________________________________________________
AliTRDtrackerFitter::AliTRDtrackerFitter(const AliTRDtrackerFitter &f)
  :TObject(f)
  ,fFitterTC(0x0)
  ,fFitterT2(0x0)
  ,fRieman1(0x0)
  ,fRieman2(0x0)
  ,fChi2TR(f.fChi2TR)
  ,fChi2TC(f.fChi2TC)
  ,fCR(f.fCR)
  ,fCC(f.fCC)
  ,fDca(f.fDca)
  ,fDzmf(f.fDzmf)
  ,fZmf(f.fZmf)
  ,fNlayers(f.fNlayers)
  ,fDebugStream(0x0)
{
  //
  // Copy Constructor (performs a deep copy) 
  //

	fRieman1  = new AliRieman(*f.fRieman1);
	fRieman2  = new AliRieman(*f.fRieman2);
	fFitterTC = new TLinearFitter(*f.fFitterTC);
	fFitterT2 = new TLinearFitter(*f.fFitterT2);
	fFitterTC->StoreData(kTRUE);
	fFitterT2->StoreData(kTRUE);
}

//_________________________________________________________________________
AliTRDtrackerFitter::~AliTRDtrackerFitter()
{
  //
  // Destructor
  //

	delete fRieman1;
	delete fRieman2;
	delete fFitterTC;
	delete fFitterT2;
}

//_________________________________________________________________________
AliTRDtrackerFitter &AliTRDtrackerFitter::operator=(const AliTRDtrackerFitter& fitter)
{
  //
  // Assignment operator using the copy Method of this class
  //
	if(this != &fitter){
		fitter.Copy(*this);
	}
	return *this;
}

//_________________________________________________________________________
void AliTRDtrackerFitter::Copy(TObject &f) const 
{
  //
  // Copy method. 
  // Performs a deep copy. Mainly used in the Assignment operator
  //

  AliTRDtrackerFitter &tf = (AliTRDtrackerFitter &) f;

	if(tf.fFitterTC)
		delete tf.fFitterTC;
	tf.fFitterTC = new TLinearFitter(*fFitterTC); 
	if(tf.fFitterT2)
		delete tf.fFitterT2;
	tf.fFitterT2 = new TLinearFitter(*fFitterT2);
	if(tf.fRieman1)
		delete tf.fRieman1;
	tf.fRieman1 = new AliRieman(*fRieman1);
	if(tf.fRieman2)
		delete fRieman2;
	tf.fRieman2 = new AliRieman(*fRieman2);
	tf.fChi2TR = fChi2TR;
	tf.fChi2TC = fChi2TC;
	tf.fCR = fCR;
	tf.fCC = fCC;
	tf.fDca =fDca;
	tf.fDzmf = fDzmf;
	tf.fZmf =fZmf;
	tf.fNlayers = fNlayers;
	tf.fDebugStream = 0x0;
	fFitterTC->StoreData(kTRUE);
	fFitterT2->StoreData(kTRUE);
}

//_________________________________________________________________________
void AliTRDtrackerFitter::FitRieman(AliTRDcluster **cl, Int_t nLayers)
{
  //
  // Performs a Rieman Fit including 4 clusters (one in each layer)
  //
	fRieman1->Reset();
	for(Int_t i = 0; i < nLayers; i++)
		fRieman1->AddPoint(cl[i]->GetX(), cl[i]->GetY(), cl[i]->GetZ(), 1, 10);
	fRieman1->Update();
}

//_________________________________________________________________________
void AliTRDtrackerFitter::FitRieman(AliTRDseedV1 *cseed, Int_t *planes)
{
  //
  // Performs a Rieman fit using four respectively six seeds
  // 2 times used: 1st with 4 parameters and Later in the full track
  // fitting Part with 6 parameters
  //

	Int_t lplanes[] = {0, 1, 2, 3, 4, 5}, *tplanes;
	Int_t nplanes;
	if(planes){
		nplanes = 4;
		tplanes = planes;
	} else {
		nplanes = 6;
		tplanes = &lplanes[0];
	}
	fRieman1->Reset();
	for(Int_t iLayer = 0; iLayer < nplanes; iLayer++){
		if(!cseed[tplanes[iLayer]].IsOK()) continue;
		fRieman1->AddPoint(cseed[tplanes[iLayer]].GetX0(), cseed[tplanes[iLayer]].GetYfitR(0), cseed[tplanes[iLayer]].GetZProb(), 1, 10);
	}
	fRieman1->Update();
}

//_________________________________________________________________________
Double_t AliTRDtrackerFitter::FitHyperplane(AliTRDseedV1 *cseed
                                          , Double_t chi2ZF
                                          , Double_t zFitter)
{
  //
  // performs a hyperplane fit
  // 
  // Two linear fitters:  one in which kz is fixed and one in which kz is a free parameter
  // checking if values are acceptable and improves the fitresults
  // Calculates curvature parameters and chisquares
  // Return: likelihood value as a measurement for the seedquality
  //

  //const Int_t kMaxTimebins = 100; // Buffer
	AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
	Int_t nTimeBins = cal->GetNumberOfTimeBins();
	Int_t npointsT = 0;
	fFitterTC->ClearPoints();
	fFitterT2->ClearPoints();
	
	Double_t fHl[6];
	Double_t xref2 = (cseed[2].GetX0() + cseed[3].GetX0())/2;
	fRieman2->Reset();
	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
		if (!cseed[iLayer].IsOK()) continue;
		fHl[iLayer] = cseed[iLayer].GetTilt();

		for (Int_t itime = 0; itime < nTimeBins; itime++) {
			if (!cseed[iLayer].IsUsable(itime)) continue;

			// X relative to the middle chamber
			Double_t x  = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
			Double_t y  = cseed[iLayer].GetY(itime);
			Double_t z  = cseed[iLayer].GetZ(itime);
			// ExB correction to the correction
			// Tilted fRieman1
			Double_t uvt[6];
			// Global x
			Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
			Double_t t  = 1.0 / (x2*x2 + y*y);
			uvt[1] = t;                 // t
			uvt[0] = 2.0 * x2 * uvt[1]; // u
			uvt[2] = 2.0 * fHl[iLayer] * uvt[1];
			uvt[3] = 2.0 * fHl[iLayer] * x * uvt[1];
			uvt[4] = 2.0 * (y + fHl[iLayer]*z) * uvt[1];
			Double_t error = 2.0 * 0.2 * uvt[1];
			fFitterT2->AddPoint(uvt,uvt[4],error);

			//Constrained fRieman1
			z = cseed[iLayer].GetZ(itime);
			uvt[0] = 2.0 * x2 * t; // u
			uvt[1] = 2.0 * fHl[iLayer] * x2 * uvt[1];
			uvt[2] = 2.0 * (y + fHl[iLayer] * (z - zFitter)) * t;	// parameter that is most propably wrong
			fFitterTC->AddPoint(uvt,uvt[2],error);
			fRieman2->AddPoint(x2,y,z,1,10);
			npointsT++;
		}
	} // Loop: iLayer

	fRieman2->Update();
	fFitterTC->Eval();
	fFitterT2->Eval();
	Double_t rpolz0 = fFitterT2->GetParameter(3);
	Double_t rpolz1 = fFitterT2->GetParameter(4);

	// Linear fitter  - not possible to make boundaries
	// Do not accept non possible z and dzdx combinations
	Bool_t acceptablez = kTRUE;
	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
		if(!cseed[iLayer].IsOK()) continue;
		Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
		if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > cseed[iLayer].GetPadLength() * 0.5 + 1.0) acceptablez = kFALSE;
	}
	if (!acceptablez) {
		fFitterT2->FixParameter(3,fRieman1->GetZat(xref2));
		fFitterT2->FixParameter(4,fRieman1->GetDZat(xref2));
		fFitterT2->Eval();
		fFitterT2->ReleaseParameter(3);
		fFitterT2->ReleaseParameter(4);
		rpolz0 = fFitterT2->GetParameter(3);
		rpolz1 = fFitterT2->GetParameter(4);
	}

	fChi2TR = fFitterT2->GetChisquare() / Float_t(npointsT);
	fChi2TC = fFitterTC->GetChisquare() / Float_t(npointsT);
	Double_t polz1c = fFitterTC->GetParameter(2);
	Double_t polz0c = polz1c * xref2;
	Double_t aC     =  fFitterTC->GetParameter(0);
	Double_t bC     =  fFitterTC->GetParameter(1);
	fCC     =  aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
	Double_t aR     =  fFitterT2->GetParameter(0);
	Double_t bR     =  fFitterT2->GetParameter(1);
	Double_t dR     =  fFitterT2->GetParameter(2);
	fCR     =  1.0 + bR*bR - dR*aR;
	fDca    =  0.0;
	if (fCR > 0.0) {
		fDca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
		fCR  =  aR / TMath::Sqrt(fCR);
	}

	
	Double_t chi2ZT2 = 0.0;
	Double_t chi2ZTC = 0.0;
	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
		if(!cseed[iLayer].IsOK()) continue;
		Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
		Double_t zTC = polz0c + polz1c * (cseed[iLayer].GetX0() - xref2);
		chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
		chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
	}
	chi2ZT2 /= TMath::Max((fNlayers - 3.0),1.0);
	chi2ZTC /= TMath::Max((fNlayers - 3.0),1.0);
	 
	
	AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
	Float_t sumdaf = 0.0;
	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
		if(!cseed[iLayer].IsOK()) continue;
		sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))/ cseed[iLayer].GetSigmaY2());
	}
	sumdaf /= Float_t (fNlayers - 2.0);

	// Cook Likelihoods
	Double_t likezf     = TMath::Exp(-chi2ZF * 0.14);
	Double_t likechi2C  = TMath::Exp(-fChi2TC * 0.677);
	Double_t likechi2TR = TMath::Exp(-fChi2TR * 0.78);
	Double_t likeaf     = TMath::Exp(-sumdaf * 3.23);
	Double_t likef = likezf * likechi2TR * likeaf;

#ifdef DEBUG
	if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 2){
		TTreeSRedirector &treeStreamer = *fDebugStream;
		treeStreamer << "FitHyperplane"
			<< "seed0.="     << &cseed[0]
			<< "seed1.="     << &cseed[1]
			<< "seed2.="     << &cseed[2]
			<< "seed3.="     << &cseed[3]
			<< "seed4.="     << &cseed[4]
			<< "seed5.="     << &cseed[5]
			<< "chi2TR="     << fChi2TR
			<< "chi2TC="     << fChi2TC
			<< "chi2ZT2="    << chi2ZT2
			<< "chi2ZTC="    << chi2ZTC
			<< "CC="         << fCC
			<< "CR="         << fCR
			<< "DR="         << dR
			<< "DCA="        << fDca
			<< "Polz0="      << polz0c
			<< "Polz1="      << polz1c
			<< "RPolz0="     << rpolz0
			<< "RPolz1="     << rpolz1
			<< "Likechi2C="  << likechi2C
			<< "Likechi2TR=" << likechi2TR
			<< "Likezf="     << likezf
			<< "Likeaf="     << likeaf
			<< "LikeF="      << likef
			<< "Rieman1.="   << fRieman1
			<< "Rieman2.="   << fRieman2
			<< "\n";
	}
#endif

	return likef;
}

//_________________________________________________________________________
void AliTRDtrackerFitter::GetHyperplaneFitChi2(Double_t *chisquares) const
{
  //
  // Getter method returns the chisquares of the hyperplane fit
  //

	chisquares[0] = fChi2TR;
	chisquares[1] = fChi2TC;
}

//_________________________________________________________________________
void AliTRDtrackerFitter::GetHyperplaneFitResults(Double_t *params) const
{
  //
  // Getter Method returning the Curvature parameters of the hyperplane fit
  //

	params[0] = fCC;
	params[1] = fCR;
	params[2] = fDca;
}

//_________________________________________________________________________
void AliTRDtrackerFitter::Reset()
{
  //
  // Resets the object
  //

	fRieman1->Reset();
	fRieman2->Reset();
	//fFitterTC->Clear();
	//fFitterT2->Clear();
}
Commit	Line	Data
e4f2f73d	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Track fitter //
	21	// //
	22	// Authors: //
	23	// Alex Bercuci <A.Bercuci@gsi.de> //
	24	// Markus Fasel <M.Fasel@gsi.de> //
	25	// //
	26	///////////////////////////////////////////////////////////////////////////////
	27
	28	#include <TROOT.h>
	29	#include <TLinearFitter.h>
	30	#include <TMath.h>
	31	#include <TTreeStream.h>
	32
	33	#include "AliLog.h"
	34
	35	#include "AliTRDseed.h"
	36	#include "AliTRDcalibDB.h"
	37	#include "AliTRDcluster.h"
	38	#include "AliTRDReconstructor.h"
	39	#include "AliTRDseedV1.h"
	40	#include "AliTRDtrackerFitter.h"
	41
	42	#define DEBUG
	43
	44	ClassImp(AliTRDtrackerFitter)
	45
	46	//_________________________________________________________________________
	47	AliTRDtrackerFitter::AliTRDtrackerFitter()
	48	:TObject()
	49	,fFitterTC(0x0)
	50	,fFitterT2(0x0)
	51	,fRieman1(0x0)
	52	,fRieman2(0x0)
	53	,fChi2TR(0)
	54	,fChi2TC(0)
	55	,fCR(0)
	56	,fCC(0)
	57	,fDca(0)
	58	,fDzmf(0)
	59	,fZmf(0)
	60	,fNlayers(0)
	61	,fDebugStream(0x0)
	62	{
	63	//
	64	// Constructor
65	//
66
67	fRieman1 = new AliRieman(1000);
68	fRieman2 = new AliRieman(1000);
69	fFitterTC = new TLinearFitter(2,"hyp2");
70	fFitterT2 = new TLinearFitter(4,"hyp4");
71	fFitterTC->StoreData(kTRUE);
72	fFitterT2->StoreData(kTRUE);
73	}
74
75	//_________________________________________________________________________
76	AliTRDtrackerFitter::AliTRDtrackerFitter(const AliTRDtrackerFitter &f)
77	:TObject(f)
78	,fFitterTC(0x0)
79	,fFitterT2(0x0)
80	,fRieman1(0x0)
81	,fRieman2(0x0)
82	,fChi2TR(f.fChi2TR)
83	,fChi2TC(f.fChi2TC)
84	,fCR(f.fCR)
85	,fCC(f.fCC)
86	,fDca(f.fDca)
87	,fDzmf(f.fDzmf)
88	,fZmf(f.fZmf)
89	,fNlayers(f.fNlayers)
90	,fDebugStream(0x0)
91	{
92	//
93	// Copy Constructor (performs a deep copy)
94	//
95
96	fRieman1 = new AliRieman(*f.fRieman1);
97	fRieman2 = new AliRieman(*f.fRieman2);
98	fFitterTC = new TLinearFitter(*f.fFitterTC);
99	fFitterT2 = new TLinearFitter(*f.fFitterT2);
100	fFitterTC->StoreData(kTRUE);
101	fFitterT2->StoreData(kTRUE);
102	}
103
104	//_________________________________________________________________________
105	AliTRDtrackerFitter::~AliTRDtrackerFitter()
106	{
107	//
108	// Destructor
109	//
110
111	delete fRieman1;
112	delete fRieman2;
113	delete fFitterTC;
114	delete fFitterT2;
115	}
116
117	//_________________________________________________________________________
118	AliTRDtrackerFitter &AliTRDtrackerFitter::operator=(const AliTRDtrackerFitter& fitter)
119	{
120	//
121	// Assignment operator using the copy Method of this class
122	//
123	if(this != &fitter){
124	fitter.Copy(*this);
125	}
126	return *this;
127	}
128
129	//_________________________________________________________________________
130	void AliTRDtrackerFitter::Copy(TObject &f) const
131	{
132	//
133	// Copy method.
134	// Performs a deep copy. Mainly used in the Assignment operator
135	//
136
137	AliTRDtrackerFitter &tf = (AliTRDtrackerFitter &) f;
138
139	if(tf.fFitterTC)
140	delete tf.fFitterTC;
141	tf.fFitterTC = new TLinearFitter(*fFitterTC);
142	if(tf.fFitterT2)
143	delete tf.fFitterT2;
144	tf.fFitterT2 = new TLinearFitter(*fFitterT2);
145	if(tf.fRieman1)
146	delete tf.fRieman1;
147	tf.fRieman1 = new AliRieman(*fRieman1);
148	if(tf.fRieman2)
149	delete fRieman2;
150	tf.fRieman2 = new AliRieman(*fRieman2);
151	tf.fChi2TR = fChi2TR;
152	tf.fChi2TC = fChi2TC;
153	tf.fCR = fCR;
154	tf.fCC = fCC;
155	tf.fDca =fDca;
156	tf.fDzmf = fDzmf;
157	tf.fZmf =fZmf;
158	tf.fNlayers = fNlayers;
159	tf.fDebugStream = 0x0;
160	fFitterTC->StoreData(kTRUE);
161	fFitterT2->StoreData(kTRUE);
162	}
163
164	//_________________________________________________________________________
165	void AliTRDtrackerFitter::FitRieman(AliTRDcluster **cl, Int_t nLayers)
166	{
167	//
168	// Performs a Rieman Fit including 4 clusters (one in each layer)
169	//
170	fRieman1->Reset();
171	for(Int_t i = 0; i < nLayers; i++)
172	fRieman1->AddPoint(cl[i]->GetX(), cl[i]->GetY(), cl[i]->GetZ(), 1, 10);
173	fRieman1->Update();
174	}
175
176	//_________________________________________________________________________
177	void AliTRDtrackerFitter::FitRieman(AliTRDseedV1 cseed, Int_t planes)
178	{
179	//
180	// Performs a Rieman fit using four respectively six seeds
181	// 2 times used: 1st with 4 parameters and Later in the full track
182	// fitting Part with 6 parameters
183	//
184
185	Int_t lplanes[] = {0, 1, 2, 3, 4, 5}, *tplanes;
186	Int_t nplanes;
187	if(planes){
188	nplanes = 4;
189	tplanes = planes;
190	} else {
191	nplanes = 6;
192	tplanes = &lplanes[0];
193	}
194	fRieman1->Reset();
195	for(Int_t iLayer = 0; iLayer < nplanes; iLayer++){
196	if(!cseed[tplanes[iLayer]].IsOK()) continue;
197	fRieman1->AddPoint(cseed[tplanes[iLayer]].GetX0(), cseed[tplanes[iLayer]].GetYfitR(0), cseed[tplanes[iLayer]].GetZProb(), 1, 10);
198	}
199	fRieman1->Update();
200	}
201
202	//_________________________________________________________________________
203	Double_t AliTRDtrackerFitter::FitHyperplane(AliTRDseedV1 *cseed
204	, Double_t chi2ZF
205	, Double_t zFitter)
206	{
207	//
208	// performs a hyperplane fit
209	//
210	// Two linear fitters: one in which kz is fixed and one in which kz is a free parameter
211	// checking if values are acceptable and improves the fitresults
212	// Calculates curvature parameters and chisquares
213	// Return: likelihood value as a measurement for the seedquality
214	//
215
216	//const Int_t kMaxTimebins = 100; // Buffer
217	AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
218	Int_t nTimeBins = cal->GetNumberOfTimeBins();
219	Int_t npointsT = 0;
220	fFitterTC->ClearPoints();
221	fFitterT2->ClearPoints();
222
223	Double_t fHl[6];
224	Double_t xref2 = (cseed[2].GetX0() + cseed[3].GetX0())/2;
225	fRieman2->Reset();
226	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
227	if (!cseed[iLayer].IsOK()) continue;
228	fHl[iLayer] = cseed[iLayer].GetTilt();
229
230	for (Int_t itime = 0; itime < nTimeBins; itime++) {
231	if (!cseed[iLayer].IsUsable(itime)) continue;
232
233	// X relative to the middle chamber
234	Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
235	Double_t y = cseed[iLayer].GetY(itime);
236	Double_t z = cseed[iLayer].GetZ(itime);
237	// ExB correction to the correction
238	// Tilted fRieman1
239	Double_t uvt[6];
240	// Global x
241	Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
242	Double_t t = 1.0 / (x2x2 + yy);
243	uvt[1] = t; // t
244	uvt[0] = 2.0 * x2 * uvt[1]; // u
245	uvt[2] = 2.0 * fHl[iLayer] * uvt[1];
246	uvt[3] = 2.0 * fHl[iLayer] * x * uvt[1];
247	uvt[4] = 2.0 * (y + fHl[iLayer]z) uvt[1];
248	Double_t error = 2.0 * 0.2 * uvt[1];
249	fFitterT2->AddPoint(uvt,uvt[4],error);
250
251	//Constrained fRieman1
252	z = cseed[iLayer].GetZ(itime);
253	uvt[0] = 2.0 * x2 * t; // u
254	uvt[1] = 2.0 * fHl[iLayer] * x2 * uvt[1];
255	uvt[2] = 2.0 * (y + fHl[iLayer] * (z - zFitter)) * t; // parameter that is most propably wrong
256	fFitterTC->AddPoint(uvt,uvt[2],error);
257	fRieman2->AddPoint(x2,y,z,1,10);
258	npointsT++;
259	}
260	} // Loop: iLayer
261
262	fRieman2->Update();
263	fFitterTC->Eval();
264	fFitterT2->Eval();
265	Double_t rpolz0 = fFitterT2->GetParameter(3);
266	Double_t rpolz1 = fFitterT2->GetParameter(4);
267
268	// Linear fitter - not possible to make boundaries
269	// Do not accept non possible z and dzdx combinations
270	Bool_t acceptablez = kTRUE;
271	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
272	if(!cseed[iLayer].IsOK()) continue;
273	Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
274	if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > cseed[iLayer].GetPadLength() * 0.5 + 1.0) acceptablez = kFALSE;
275	}
276	if (!acceptablez) {
277	fFitterT2->FixParameter(3,fRieman1->GetZat(xref2));
278	fFitterT2->FixParameter(4,fRieman1->GetDZat(xref2));
279	fFitterT2->Eval();
280	fFitterT2->ReleaseParameter(3);
281	fFitterT2->ReleaseParameter(4);
282	rpolz0 = fFitterT2->GetParameter(3);
283	rpolz1 = fFitterT2->GetParameter(4);
284	}
285
286	fChi2TR = fFitterT2->GetChisquare() / Float_t(npointsT);
287	fChi2TC = fFitterTC->GetChisquare() / Float_t(npointsT);
288	Double_t polz1c = fFitterTC->GetParameter(2);
289	Double_t polz0c = polz1c * xref2;
290	Double_t aC = fFitterTC->GetParameter(0);
291	Double_t bC = fFitterTC->GetParameter(1);
292	fCC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
293	Double_t aR = fFitterT2->GetParameter(0);
294	Double_t bR = fFitterT2->GetParameter(1);
295	Double_t dR = fFitterT2->GetParameter(2);
296	fCR = 1.0 + bRbR - dRaR;
297	fDca = 0.0;
298	if (fCR > 0.0) {
299	fDca = -dR / (TMath::Sqrt(1.0 + bRbR - dRaR) + TMath::Sqrt(1.0 + bR*bR));
300	fCR = aR / TMath::Sqrt(fCR);
301	}
302
303
304	Double_t chi2ZT2 = 0.0;
305	Double_t chi2ZTC = 0.0;
306	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
307	if(!cseed[iLayer].IsOK()) continue;
308	Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
309	Double_t zTC = polz0c + polz1c * (cseed[iLayer].GetX0() - xref2);
310	chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
311	chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
312	}
313	chi2ZT2 /= TMath::Max((fNlayers - 3.0),1.0);
314	chi2ZTC /= TMath::Max((fNlayers - 3.0),1.0);
315
316
317	AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
318	Float_t sumdaf = 0.0;
319	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
320	if(!cseed[iLayer].IsOK()) continue;
321	sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))/ cseed[iLayer].GetSigmaY2());
322	}
323	sumdaf /= Float_t (fNlayers - 2.0);
324
325	// Cook Likelihoods
326	Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
327	Double_t likechi2C = TMath::Exp(-fChi2TC * 0.677);
328	Double_t likechi2TR = TMath::Exp(-fChi2TR * 0.78);
329	Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
330	Double_t likef = likezf * likechi2TR * likeaf;
331
332	#ifdef DEBUG
333	if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 2){
334	TTreeSRedirector &treeStreamer = *fDebugStream;
335	treeStreamer << "FitHyperplane"
336	<< "seed0.=" << &cseed[0]
337	<< "seed1.=" << &cseed[1]
338	<< "seed2.=" << &cseed[2]
339	<< "seed3.=" << &cseed[3]
340	<< "seed4.=" << &cseed[4]
341	<< "seed5.=" << &cseed[5]
342	<< "chi2TR=" << fChi2TR
343	<< "chi2TC=" << fChi2TC
344	<< "chi2ZT2=" << chi2ZT2
345	<< "chi2ZTC=" << chi2ZTC
346	<< "CC=" << fCC
347	<< "CR=" << fCR
348	<< "DR=" << dR
349	<< "DCA=" << fDca
350	<< "Polz0=" << polz0c
351	<< "Polz1=" << polz1c
352	<< "RPolz0=" << rpolz0
353	<< "RPolz1=" << rpolz1
354	<< "Likechi2C=" << likechi2C
355	<< "Likechi2TR=" << likechi2TR
356	<< "Likezf=" << likezf
357	<< "Likeaf=" << likeaf
358	<< "LikeF=" << likef
359	<< "Rieman1.=" << fRieman1
360	<< "Rieman2.=" << fRieman2
361	<< "\n";
362	}
363	#endif
364
365	return likef;
366	}
367
368	//_________________________________________________________________________
369	void AliTRDtrackerFitter::GetHyperplaneFitChi2(Double_t *chisquares) const
370	{
371	//
372	// Getter method returns the chisquares of the hyperplane fit
373	//
374
375	chisquares[0] = fChi2TR;
376	chisquares[1] = fChi2TC;
377	}
378
379	//_________________________________________________________________________
380	void AliTRDtrackerFitter::GetHyperplaneFitResults(Double_t *params) const
381	{
382	//
383	// Getter Method returning the Curvature parameters of the hyperplane fit
384	//
385
386	params[0] = fCC;
387	params[1] = fCR;
388	params[2] = fDca;
389	}
390
391	//_________________________________________________________________________
392	void AliTRDtrackerFitter::Reset()
393	{
394	//
395	// Resets the object
396	//
397
398	fRieman1->Reset();
399	fRieman2->Reset();
400	//fFitterTC->Clear();
401	//fFitterT2->Clear();
402	}