[u/mrichter/AliRoot.git] / TRD / AliTRDgtuTrack.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.                  *
 **************************************************************************/

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//                                                                           //
//  TRD module global track (GTU)                                            //
//                                                                           //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TMatrixD.h>
#include <TObjArray.h>

#include "AliTRDgeometry.h"
#include "AliTRDcalibDB.h"
#include "Cal/AliTRDCalPIDLQ.h"

#include "AliTRDltuTracklet.h"

#include "AliTRDgtuTrack.h"

ClassImp(AliTRDgtuTrack)

//_____________________________________________________________________________
AliTRDgtuTrack::AliTRDgtuTrack()
{
  //
  // Default constructor
  //

  fYproj = 0.0;
  fZproj = 0.0;
  fSlope = 0.0;
  fDetector = -1;
  fNtracklets = 0;
  fNplanes = 0;
  fNclusters = 0;
  fPt = 0.0;
  fPhi = 0.0;
  fEta = 0.0;
  fLabel = -1;
  fPID = 0.0;
  fIsElectron = kFALSE;

  fTracklets = new TObjArray(400);

}

//_____________________________________________________________________________
AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack& track):
  TObject(track),
  fTracklets(NULL),
  fYproj(track.fYproj),
  fZproj(track.fZproj),
  fSlope(track.fSlope),
  fDetector(track.fDetector),
  fNtracklets(track.fNtracklets),
  fNplanes(track.fNplanes),
  fNclusters(track.fNclusters),
  fPt(track.fPt),
  fPhi(track.fPhi),
  fEta(track.fEta),
  fLabel(track.fLabel),
  fPID(track.fPID),
  fIsElectron(track.fIsElectron)
{
  //
  // copy contructor
  //

}

//_____________________________________________________________________________
AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
{
  //
  // assignment operator
  //

  if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this); 
  return *this;

}

//_____________________________________________________________________________
void AliTRDgtuTrack::Copy(TObject &t) const
{
  //
  // copy function
  //

  ((AliTRDgtuTrack &) t).fTracklets  = NULL;
  ((AliTRDgtuTrack &) t).fYproj      = fYproj;
  ((AliTRDgtuTrack &) t).fZproj      = fZproj;
  ((AliTRDgtuTrack &) t).fSlope      = fSlope;
  ((AliTRDgtuTrack &) t).fDetector   = fDetector;
  ((AliTRDgtuTrack &) t).fNtracklets = fNtracklets;
  ((AliTRDgtuTrack &) t).fNplanes    = fNplanes;
  ((AliTRDgtuTrack &) t).fNclusters  = fNclusters;
  ((AliTRDgtuTrack &) t).fPt         = fPt;
  ((AliTRDgtuTrack &) t).fPhi        = fPhi;
  ((AliTRDgtuTrack &) t).fEta        = fEta;
  ((AliTRDgtuTrack &) t).fLabel      = fLabel;
  ((AliTRDgtuTrack &) t).fPID        = fPID;
  ((AliTRDgtuTrack &) t).fIsElectron = fIsElectron;

}

//_____________________________________________________________________________
AliTRDgtuTrack::~AliTRDgtuTrack()
{
  //
  // destructor
  //

}

//_____________________________________________________________________________
void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
{
  //
  // Add a LTU tracklet to this track
  //

  Tracklets()->Add(trk);

}

//_____________________________________________________________________________
AliTRDltuTracklet* AliTRDgtuTrack::GetTracklet(Int_t pos) const
{
  //
  // Return LTU tracklet at position "pos"
  //

  if (fTracklets == 0) return 0;
  void * trk = fTracklets->UncheckedAt(pos);
  if (trk == 0) return 0;

  return (AliTRDltuTracklet*)trk;

}

//_____________________________________________________________________________
Int_t AliTRDgtuTrack::Compare(const TObject * o) const
{
  //
  // Compare function for sorting the tracks
  //

  AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack*)o;

  if (fYproj <  gtutrack->GetYproj()) return -1;
  if (fYproj == gtutrack->GetYproj()) return  0;

  return +1;

}

//_____________________________________________________________________________
void AliTRDgtuTrack::Reset()
{
  //
  // Reset the track information
  //

  fYproj = 0.0;
  fZproj = 0.0;
  fSlope = 0.0;
  fDetector = -1;
  fNtracklets = 0;
  fNplanes = 0;
  fNclusters = 0;
  fPt = 0.0;
  fPhi = 0.0;
  fEta = 0.0;
  fLabel = -1;
  fPID = 0.0;
  fIsElectron = kFALSE;
  
}

//_____________________________________________________________________________
void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
{
  //
  // Calculate the kinematics of the found track
  //

  AliTRDltuTracklet *trk;
  Int_t nTracklets = GetNtracklets();
  Float_t fC[kNmaxTrk][3];            // X, Y, Z  coordinates of segments

  fYproj = 0.0;
  fZproj = 0.0;
  fSlope = 0.0;
  fNclusters = 0;
  fNplanes = 0;
  fNtracklets = GetNtracklets();
  Int_t inDetector[kNplan];
  for (Int_t i = 0; i < kNplan; i++) inDetector[i] = -1;
  Int_t iDet, nDet = 0;
  Bool_t newDetector;
  for (Int_t i = 0; i < nTracklets; i++) {

    trk = GetTracklet(i);
    fYproj += trk->GetYproj(xpl);
    fZproj += trk->GetZproj(xpl);
    fSlope += trk->GetSlope();
    fNclusters += trk->GetNclusters();
    iDet = trk->GetDetector();

    newDetector = kTRUE;
    for (Int_t id = 0; id < nDet; id++) {
      if (iDet == inDetector[id]) {
	newDetector = kFALSE;
	break;
      }
    }
    if (newDetector) {
      inDetector[nDet++] = iDet;
      fNplanes++;
    }

    fC[i][0] = trk->GetTime0();
    fC[i][1] = trk->GetOffset();
    fC[i][2] = trk->GetRowz();

  }
  fYproj /= (Float_t)nTracklets;
  fZproj /= (Float_t)nTracklets;
  fSlope /= (Float_t)nTracklets;

  Float_t x[kNmaxTrk+1], y[kNmaxTrk+1], z[kNmaxTrk+1];
  Bool_t count[kNmaxTrk];
  for (Int_t i = 0; i < kNmaxTrk; i++) count[i] = kFALSE;

  Int_t iXmin = -1;
  Int_t j = 0;
  x[0] = y[0] = z[0] = 0.0;
  while (j < nTracklets) {
    iXmin = -1;
    for (Int_t i = 0; i < nTracklets; i++) {
      if (count[i]) continue;
      if (iXmin == -1) {
	iXmin = i;
	continue;
      }
      if (fC[i][0] < fC[iXmin][0]) iXmin = i;
    }
    x[j+1] = fC[iXmin][0];
    y[j+1] = fC[iXmin][1];
    z[j+1] = fC[iXmin][2];
    j++;
    count[iXmin] = kTRUE;
  }

  TMatrixD smatrix(2,2);
  TMatrixD sums(2,1);
  TMatrixD res(2,1);
  Double_t xv, yv;
  Float_t a, b;

  smatrix.Zero();
  sums.Zero();
  for (Int_t i = 0; i < nTracklets; i++) {
    xv = (Double_t)x[i+1];
    yv = (Double_t)y[i+1];
    smatrix(0,0) += 1.0;
    smatrix(1,1) += xv*xv;
    smatrix(0,1) += xv;
    smatrix(1,0) += xv;
    sums(0,0)    += yv;
    sums(1,0)    += xv*yv;
  }
  res = smatrix.Invert() * sums;
  a = res(0,0);
  b = res(1,0);

  Float_t dist = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);

  Float_t fx1 = x[1]          + dist * (Float_t)(nTracklets-1)/6.0;
  Float_t fy1 = a + b * fx1;
  Float_t fx2 = x[nTracklets] - dist * (Float_t)(nTracklets-1)/6.0;
  Float_t fy2 = a + b * fx2;
  Float_t d12 = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
  Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
  Float_t r = (d12/2.0)/TMath::Sin(alpha);

  fPt = 0.3 * field * 0.01 * r;

  Float_t d1 = fx1*fx1+fy1*fy1;
  Float_t d2 = fx2*fx2+fy2*fy2;
  Float_t d  = fx1*fy2-fx2*fy1;
  
  Float_t xc = (d1*fy2-d2*fy1)/(2*d);
  Float_t yc = (d2*fx1-d1*fx2)/(2*d);

  if (yc != 0.0) {
    fPhi = TMath::ATan(xc/yc);
  } else {
    fPhi = TMath::PiOver2();
  }

  fPhi *= 180.0/TMath::Pi();

  smatrix.Zero();
  sums.Zero();
  for (Int_t i = 0; i < nTracklets+1; i++) {
    xv = (Double_t)z[i];
    yv = (Double_t)x[i];
    smatrix(0,0) += 1.0;
    smatrix(1,1) += xv*xv;
    smatrix(0,1) += xv;
    smatrix(1,0) += xv;
    sums(0,0)    += yv;
    sums(1,0)    += xv*yv;
  }
  res = smatrix.Invert() * sums;
  a = res(0,0);
  b = res(1,0);
  Float_t theta = TMath::ATan(b);
  
  if (theta < 0.0) theta = TMath::Pi() + theta;
  
  if (theta == 0.0) {
    fEta = 0.0;
  } else {
    fEta = -TMath::Log(TMath::Tan(theta/2.0));
  }
  
}

//_____________________________________________________________________________
void AliTRDgtuTrack::MakePID()
{
  //
  // Electron likelihood signal
  //

  AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
  if (!calibration)
  {
    Error("MakePID","No instance of AliTRDcalibDB.");
    return;  
  }
  const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
  
  AliTRDltuTracklet *trk;
  Int_t nTracklets = GetNtracklets();
  Int_t det, pla;
  Float_t sl, th, q, probPio = 1.0, probEle = 1.0;
  for (Int_t i = 0; i < nTracklets; i++) {

    trk = GetTracklet(i);

    sl = TMath::Abs(trk->GetSlope());     // tracklet inclination in X-y plane
    th = trk->GetRowz()/trk->GetTime0();  // tracklet inclination in X-z plane
    th = TMath::ATan(TMath::Abs(th));

    q = trk->GetQ() 
      * TMath::Cos(sl/180.0*TMath::Pi()) 
      * TMath::Cos(th/180.0*TMath::Pi());

    det = trk->GetDetector();
    pla = trk->GetPlane(det);

    // unclear yet factor to match the PS distributions = 5.8
    // not explained only by the tail filter ...

    // AliRoot v4-03-07 , v4-03-Release
    //q = q * 5.8;

    // HEAD28Mar06
    // Temporary (B. Vulpescu):
    // The charge distributions do not match the new changes in simulation (A. Bercuci),
    // which are nevertheless now in agreement with the beam tests.
    // Some tricks will be used to still have reasonable results
    // To match the existing charge distributions, the charge per layer has to be modified
    // as folows:
    /*
      if (k == 0) {
      // electrons
      q = 4.3 * q + 95.0;
      } else {
      // others
      q = 4.2 * q + 70.0;
      }
    */
    // Since at tracking time we have no information on the particle type, we will modify
    // instead the charge distributions accordingly. This will slow down the sampling.
    // The modified distributions are in TRDdEdxHistogramsV1_BV.root and the CDB has
    // been regenerated with AliTRDCreateDummyCDB.C
    // The new PIDLQ data base has the version :
    // I-AliCDBLocal::Get: CDB object retrieved: path "TRD/Calib/PIDLQ"; run range [0,0];
    // version v0_s1

    // HEAD28Mar06 new distributions: factor = 6.0

    probEle *= pd->GetProbability(0,TMath::Abs(fPt),q*6.0);
    probPio *= pd->GetProbability(2,TMath::Abs(fPt),q*6.0);

  }

  if ((probEle+probPio) > 0.0) {
    fPID = probEle/(probEle+probPio);
  } else {
    fPID = 0.0;
  }

  // Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
  // P [GeV/c]  fPIDcut (between 0 and 1)
  // 2          0.925
  // 3          0.915
  // 4          0.875
  // 5          0.855
  // 6          0.845
  // 8          0.785
  // 10         0.735
  //
  // PIDcut = 0.978 - 0.024 * P[GeV/c]
  //Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);

  // HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
  //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);

  // HEAD28Mar06 with new generated distributions (pol2 fit)
  Float_t absPt = TMath::Abs(fPt);
  //Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt;  // 800 bins in dEdx
  Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt;    // 200 bins in dEdx

  fIsElectron = kFALSE;
  if (fPID >= fPIDcut) {
    fIsElectron = kTRUE;
  }

}

//_____________________________________________________________________________
void AliTRDgtuTrack::CookLabel()
{
  //
  // Cook the track label from tracklets labels
  //

  AliTRDltuTracklet *trk;

  const Int_t kMaxTracks = 10;
  Int_t trackLabel[kMaxTracks];
  Int_t trackCount[kMaxTracks];
  for (Int_t it = 0; it < kMaxTracks; it++) {
    trackLabel[it] = -1;
    trackCount[it] = 0;
  }

  Bool_t counted;
  Int_t label, nTracks = 0;
  for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {

    trk = GetTracklet(itrk);

    if (trk->GetLabel() == -1) continue;

    label = trk->GetLabel();

    counted = kFALSE;
    for (Int_t it = 0; it < nTracks; it++) {
      if (label == trackLabel[it]) {
	trackCount[it]++;
	counted = kTRUE;
	break;
      }
    }
    if (!counted) {
      trackLabel[nTracks] = label;
      trackCount[nTracks]++;
      nTracks++;
      if (nTracks == kMaxTracks) {
	Warning("CookLabel","Too many tracks for this tracklet.");
	nTracks--;
	break;
      }
    }
    
  }

  Float_t frac = 4.0/5.0;
  for (Int_t it = 0; it < kMaxTracks; it++) {
    if (trackCount[it] >= (Int_t)(frac*fNtracklets)) {
      fLabel = trackLabel[it];
      break;
    }
  }

}
Commit	Line	Data
e3b2b5e5	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	///////////////////////////////////////////////////////////////////////////////
	17	// //
	18	// //
	19	// TRD module global track (GTU) //
	20	// //
	21	// //
	22	///////////////////////////////////////////////////////////////////////////////
	23
	24	#include <TMatrixD.h>
	25	#include <TObjArray.h>
	26
	27	#include "AliTRDgeometry.h"
	28	#include "AliTRDcalibDB.h"
	29	#include "Cal/AliTRDCalPIDLQ.h"
	30
	31	#include "AliTRDltuTracklet.h"
	32
	33	#include "AliTRDgtuTrack.h"
	34
	35	ClassImp(AliTRDgtuTrack)
	36
	37	//_____________________________________________________________________________
	38	AliTRDgtuTrack::AliTRDgtuTrack()
	39	{
	40	//
	41	// Default constructor
	42	//
	43
	44	fYproj = 0.0;
	45	fZproj = 0.0;
	46	fSlope = 0.0;
	47	fDetector = -1;
	48	fNtracklets = 0;
	49	fNplanes = 0;
	50	fNclusters = 0;
	51	fPt = 0.0;
	52	fPhi = 0.0;
	53	fEta = 0.0;
	54	fLabel = -1;
	55	fPID = 0.0;
	56	fIsElectron = kFALSE;
	57
	58	fTracklets = new TObjArray(400);
	59
	60	}
	61
	62	//_____________________________________________________________________________
	63	AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack& track):
	64	TObject(track),
65	fTracklets(NULL),
66	fYproj(track.fYproj),
67	fZproj(track.fZproj),
68	fSlope(track.fSlope),
69	fDetector(track.fDetector),
70	fNtracklets(track.fNtracklets),
71	fNplanes(track.fNplanes),
72	fNclusters(track.fNclusters),
73	fPt(track.fPt),
74	fPhi(track.fPhi),
75	fEta(track.fEta),
76	fLabel(track.fLabel),
77	fPID(track.fPID),
78	fIsElectron(track.fIsElectron)
79	{
80	//
81	// copy contructor
82	//
83
84	}
85
86	//_____________________________________________________________________________
87	AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
88	{
89	//
90	// assignment operator
91	//
92
93	if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this);
94	return *this;
95
96	}
97
98	//_____________________________________________________________________________
99	void AliTRDgtuTrack::Copy(TObject &t) const
100	{
101	//
102	// copy function
103	//
104
105	((AliTRDgtuTrack &) t).fTracklets = NULL;
106	((AliTRDgtuTrack &) t).fYproj = fYproj;
107	((AliTRDgtuTrack &) t).fZproj = fZproj;
108	((AliTRDgtuTrack &) t).fSlope = fSlope;
109	((AliTRDgtuTrack &) t).fDetector = fDetector;
110	((AliTRDgtuTrack &) t).fNtracklets = fNtracklets;
111	((AliTRDgtuTrack &) t).fNplanes = fNplanes;
112	((AliTRDgtuTrack &) t).fNclusters = fNclusters;
113	((AliTRDgtuTrack &) t).fPt = fPt;
114	((AliTRDgtuTrack &) t).fPhi = fPhi;
115	((AliTRDgtuTrack &) t).fEta = fEta;
116	((AliTRDgtuTrack &) t).fLabel = fLabel;
117	((AliTRDgtuTrack &) t).fPID = fPID;
118	((AliTRDgtuTrack &) t).fIsElectron = fIsElectron;
119
120	}
121
122	//_____________________________________________________________________________
123	AliTRDgtuTrack::~AliTRDgtuTrack()
124	{
125	//
126	// destructor
127	//
128
129	}
130
131	//_____________________________________________________________________________
132	void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
133	{
134	//
135	// Add a LTU tracklet to this track
136	//
137
138	Tracklets()->Add(trk);
139
140	}
141
142	//_____________________________________________________________________________
143	AliTRDltuTracklet* AliTRDgtuTrack::GetTracklet(Int_t pos) const
144	{
145	//
146	// Return LTU tracklet at position "pos"
147	//
148
149	if (fTracklets == 0) return 0;
150	void * trk = fTracklets->UncheckedAt(pos);
151	if (trk == 0) return 0;
152
153	return (AliTRDltuTracklet*)trk;
154
155	}
156
157	//_____________________________________________________________________________
158	Int_t AliTRDgtuTrack::Compare(const TObject * o) const
159	{
160	//
161	// Compare function for sorting the tracks
162	//
163
164	AliTRDgtuTrack gtutrack = (AliTRDgtuTrack)o;
165
166	if (fYproj < gtutrack->GetYproj()) return -1;
167	if (fYproj == gtutrack->GetYproj()) return 0;
168
169	return +1;
170
171	}
172
173	//_____________________________________________________________________________
174	void AliTRDgtuTrack::Reset()
175	{
176	//
177	// Reset the track information
178	//
179
180	fYproj = 0.0;
181	fZproj = 0.0;
182	fSlope = 0.0;
183	fDetector = -1;
184	fNtracklets = 0;
185	fNplanes = 0;
186	fNclusters = 0;
187	fPt = 0.0;
188	fPhi = 0.0;
189	fEta = 0.0;
190	fLabel = -1;
191	fPID = 0.0;
192	fIsElectron = kFALSE;
193
194	}
195
196	//_____________________________________________________________________________
197	void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
198	{
199	//
200	// Calculate the kinematics of the found track
201	//
202
203	AliTRDltuTracklet *trk;
204	Int_t nTracklets = GetNtracklets();
205	Float_t fC[kNmaxTrk][3]; // X, Y, Z coordinates of segments
206
207	fYproj = 0.0;
208	fZproj = 0.0;
209	fSlope = 0.0;
210	fNclusters = 0;
211	fNplanes = 0;
212	fNtracklets = GetNtracklets();
213	Int_t inDetector[kNplan];
214	for (Int_t i = 0; i < kNplan; i++) inDetector[i] = -1;
215	Int_t iDet, nDet = 0;
216	Bool_t newDetector;
217	for (Int_t i = 0; i < nTracklets; i++) {
218
219	trk = GetTracklet(i);
220	fYproj += trk->GetYproj(xpl);
221	fZproj += trk->GetZproj(xpl);
222	fSlope += trk->GetSlope();
223	fNclusters += trk->GetNclusters();
224	iDet = trk->GetDetector();
225
226	newDetector = kTRUE;
227	for (Int_t id = 0; id < nDet; id++) {
228	if (iDet == inDetector[id]) {
229	newDetector = kFALSE;
230	break;
231	}
232	}
233	if (newDetector) {
234	inDetector[nDet++] = iDet;
235	fNplanes++;
236	}
237
238	fC[i][0] = trk->GetTime0();
239	fC[i][1] = trk->GetOffset();
240	fC[i][2] = trk->GetRowz();
241
242	}
243	fYproj /= (Float_t)nTracklets;
244	fZproj /= (Float_t)nTracklets;
245	fSlope /= (Float_t)nTracklets;
246
247	Float_t x[kNmaxTrk+1], y[kNmaxTrk+1], z[kNmaxTrk+1];
248	Bool_t count[kNmaxTrk];
249	for (Int_t i = 0; i < kNmaxTrk; i++) count[i] = kFALSE;
250
251	Int_t iXmin = -1;
252	Int_t j = 0;
253	x[0] = y[0] = z[0] = 0.0;
254	while (j < nTracklets) {
255	iXmin = -1;
256	for (Int_t i = 0; i < nTracklets; i++) {
257	if (count[i]) continue;
258	if (iXmin == -1) {
259	iXmin = i;
260	continue;
261	}
262	if (fC[i][0] < fC[iXmin][0]) iXmin = i;
263	}
264	x[j+1] = fC[iXmin][0];
265	y[j+1] = fC[iXmin][1];
266	z[j+1] = fC[iXmin][2];
267	j++;
268	count[iXmin] = kTRUE;
269	}
270
271	TMatrixD smatrix(2,2);
272	TMatrixD sums(2,1);
273	TMatrixD res(2,1);
274	Double_t xv, yv;
275	Float_t a, b;
276
277	smatrix.Zero();
278	sums.Zero();
279	for (Int_t i = 0; i < nTracklets; i++) {
280	xv = (Double_t)x[i+1];
281	yv = (Double_t)y[i+1];
282	smatrix(0,0) += 1.0;
283	smatrix(1,1) += xv*xv;
284	smatrix(0,1) += xv;
285	smatrix(1,0) += xv;
286	sums(0,0) += yv;
287	sums(1,0) += xv*yv;
288	}
289	res = smatrix.Invert() * sums;
290	a = res(0,0);
291	b = res(1,0);
292
293	Float_t dist = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
294
295	Float_t fx1 = x[1] + dist * (Float_t)(nTracklets-1)/6.0;
296	Float_t fy1 = a + b * fx1;
297	Float_t fx2 = x[nTracklets] - dist * (Float_t)(nTracklets-1)/6.0;
298	Float_t fy2 = a + b * fx2;
299	Float_t d12 = TMath::Sqrt((fx2-fx1)(fx2-fx1)+(fy2-fy1)(fy2-fy1));
300	Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
301	Float_t r = (d12/2.0)/TMath::Sin(alpha);
302
303	fPt = 0.3 * field * 0.01 * r;
304
305	Float_t d1 = fx1fx1+fy1fy1;
306	Float_t d2 = fx2fx2+fy2fy2;
307	Float_t d = fx1fy2-fx2fy1;
308
309	Float_t xc = (d1fy2-d2fy1)/(2*d);
310	Float_t yc = (d2fx1-d1fx2)/(2*d);
311
312	if (yc != 0.0) {
313	fPhi = TMath::ATan(xc/yc);
314	} else {
315	fPhi = TMath::PiOver2();
316	}
317
318	fPhi *= 180.0/TMath::Pi();
319
320	smatrix.Zero();
321	sums.Zero();
322	for (Int_t i = 0; i < nTracklets+1; i++) {
323	xv = (Double_t)z[i];
324	yv = (Double_t)x[i];
325	smatrix(0,0) += 1.0;
326	smatrix(1,1) += xv*xv;
327	smatrix(0,1) += xv;
328	smatrix(1,0) += xv;
329	sums(0,0) += yv;
330	sums(1,0) += xv*yv;
331	}
332	res = smatrix.Invert() * sums;
333	a = res(0,0);
334	b = res(1,0);
335	Float_t theta = TMath::ATan(b);
336
337	if (theta < 0.0) theta = TMath::Pi() + theta;
338
339	if (theta == 0.0) {
340	fEta = 0.0;
341	} else {
342	fEta = -TMath::Log(TMath::Tan(theta/2.0));
343	}
344
345	}
346
347	//_____________________________________________________________________________
348	void AliTRDgtuTrack::MakePID()
349	{
350	//
351	// Electron likelihood signal
352	//
353
354	AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
355	if (!calibration)
356	{
357	Error("MakePID","No instance of AliTRDcalibDB.");
358	return;
359	}
360	const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
361
362	AliTRDltuTracklet *trk;
363	Int_t nTracklets = GetNtracklets();
364	Int_t det, pla;
365	Float_t sl, th, q, probPio = 1.0, probEle = 1.0;
366	for (Int_t i = 0; i < nTracklets; i++) {
367
368	trk = GetTracklet(i);
369
370	sl = TMath::Abs(trk->GetSlope()); // tracklet inclination in X-y plane
371	th = trk->GetRowz()/trk->GetTime0(); // tracklet inclination in X-z plane
372	th = TMath::ATan(TMath::Abs(th));
373
374	q = trk->GetQ()
375	* TMath::Cos(sl/180.0*TMath::Pi())
376	* TMath::Cos(th/180.0*TMath::Pi());
377
378	det = trk->GetDetector();
379	pla = trk->GetPlane(det);
380
381	// unclear yet factor to match the PS distributions = 5.8
382	// not explained only by the tail filter ...
383
384	// AliRoot v4-03-07 , v4-03-Release
385	//q = q * 5.8;
386
4c8471bd	387	// HEAD28Mar06
e3b2b5e5	388	// Temporary (B. Vulpescu):
	389	// The charge distributions do not match the new changes in simulation (A. Bercuci),
	390	// which are nevertheless now in agreement with the beam tests.
	391	// Some tricks will be used to still have reasonable results
	392	// To match the existing charge distributions, the charge per layer has to be modified
	393	// as folows:
	394	/*
	395	if (k == 0) {
	396	// electrons
	397	q = 4.3 * q + 95.0;
	398	} else {
	399	// others
	400	q = 4.2 * q + 70.0;
	401	}
	402	*/
	403	// Since at tracking time we have no information on the particle type, we will modify
	404	// instead the charge distributions accordingly. This will slow down the sampling.
	405	// The modified distributions are in TRDdEdxHistogramsV1_BV.root and the CDB has
	406	// been regenerated with AliTRDCreateDummyCDB.C
	407	// The new PIDLQ data base has the version :
	408	// I-AliCDBLocal::Get: CDB object retrieved: path "TRD/Calib/PIDLQ"; run range [0,0];
	409	// version v0_s1
	410
4c8471bd	411	// HEAD28Mar06 new distributions: factor = 6.0
	412
	413	probEle = pd->GetProbability(0,TMath::Abs(fPt),q6.0);
	414	probPio = pd->GetProbability(2,TMath::Abs(fPt),q6.0);
e3b2b5e5	415
	416	}
	417
	418	if ((probEle+probPio) > 0.0) {
	419	fPID = probEle/(probEle+probPio);
	420	} else {
	421	fPID = 0.0;
	422	}
	423
	424	// Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
	425	// P [GeV/c] fPIDcut (between 0 and 1)
	426	// 2 0.925
	427	// 3 0.915
	428	// 4 0.875
	429	// 5 0.855
	430	// 6 0.845
	431	// 8 0.785
	432	// 10 0.735
	433	//
	434	// PIDcut = 0.978 - 0.024 * P[GeV/c]
	435	//Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);
	436
	437	// HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
4c8471bd	438	//Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
	439
	440	// HEAD28Mar06 with new generated distributions (pol2 fit)
	441	Float_t absPt = TMath::Abs(fPt);
	442	//Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt; // 800 bins in dEdx
	443	Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt; // 200 bins in dEdx
e3b2b5e5	444
	445	fIsElectron = kFALSE;
	446	if (fPID >= fPIDcut) {
	447	fIsElectron = kTRUE;
	448	}
	449
	450	}
	451
	452	//_____________________________________________________________________________
	453	void AliTRDgtuTrack::CookLabel()
	454	{
	455	//
	456	// Cook the track label from tracklets labels
	457	//
	458
	459	AliTRDltuTracklet *trk;
	460
	461	const Int_t kMaxTracks = 10;
	462	Int_t trackLabel[kMaxTracks];
	463	Int_t trackCount[kMaxTracks];
	464	for (Int_t it = 0; it < kMaxTracks; it++) {
	465	trackLabel[it] = -1;
	466	trackCount[it] = 0;
	467	}
	468
	469	Bool_t counted;
	470	Int_t label, nTracks = 0;
	471	for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
	472
	473	trk = GetTracklet(itrk);
	474
	475	if (trk->GetLabel() == -1) continue;
	476
	477	label = trk->GetLabel();
	478
	479	counted = kFALSE;
	480	for (Int_t it = 0; it < nTracks; it++) {
	481	if (label == trackLabel[it]) {
	482	trackCount[it]++;
	483	counted = kTRUE;
	484	break;
	485	}
	486	}
	487	if (!counted) {
	488	trackLabel[nTracks] = label;
	489	trackCount[nTracks]++;
	490	nTracks++;
	491	if (nTracks == kMaxTracks) {
	492	Warning("CookLabel","Too many tracks for this tracklet.");
	493	nTracks--;
	494	break;
	495	}
	496	}
	497
	498	}
	499
	500	Float_t frac = 4.0/5.0;
	501	for (Int_t it = 0; it < kMaxTracks; it++) {
	502	if (trackCount[it] >= (Int_t)(frac*fNtracklets)) {
	503	fLabel = trackLabel[it];
	504	break;
	505	}
	506	}
	507
508	}
509