1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 #include "AliTracker.h"
22 #include "AliTRDgeometry.h"
23 #include "AliTRDcluster.h"
24 #include "AliTRDtrack.h"
25 #include "AliTRDtracklet.h"
26 #include "AliTRDcalibDB.h"
30 ///////////////////////////////////////////////////////////////////////////////
32 // Represents a reconstructed TRD track //
33 // Local TRD Kalman track //
34 // Part of the old TRD tracking code //
36 ///////////////////////////////////////////////////////////////////////////////
38 //_____________________________________________________________________________
39 AliTRDtrack::AliTRDtrack()
45 ,fClusterOwner(kFALSE)
60 // AliTRDtrack default constructor
63 for (Int_t i = 0; i < kNplane; i++) {
64 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
65 fdEdxPlane[i][j] = 0.0;
71 fTrackletIndex[i] = -1;
73 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
74 fPID[ispec] = 1.0 / AliPID::kSPECIES;
77 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
84 for (Int_t i = 0; i < 3; i++) {
90 //_____________________________________________________________________________
91 AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index
92 , const Double_t p[5], const Double_t cov[15]
93 , Double_t x, Double_t alpha)
99 ,fClusterOwner(kFALSE)
114 // The main AliTRDtrack constructor.
117 Double_t cnv = 1.0 / (GetBz() * kB2C);
119 Double_t pp[5] = { p[0]
125 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
126 Double_t c32 = x*cov[13] - cov[ 8];
127 Double_t c20 = x*cov[10] - cov[ 3];
128 Double_t c21 = x*cov[11] - cov[ 4];
129 Double_t c42 = x*cov[14] - cov[12];
131 Double_t cc[15] = { cov[ 0]
134 , cov[ 6], cov[ 7], c32, cov[ 9]
135 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
138 SetNumberOfClusters(1);
142 for (Int_t i = 0; i < kNplane; i++) {
143 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
144 fdEdxPlane[i][j] = 0.0;
146 fTimBinPlane[i] = -1;
150 fTrackletIndex[i] = -1;
152 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
153 fPID[ispec] = 1.0 / AliPID::kSPECIES;
156 Double_t q = TMath::Abs(c->GetQ());
157 Double_t s = GetSnp();
158 Double_t t = GetTgl();
160 q *= TMath::Sqrt((1-s*s)/(1+t*t));
164 for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
171 for (Int_t i = 0; i < 3;i++) {
177 //_____________________________________________________________________________
178 AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
180 ,fSeedLab(t.GetSeedLabel())
183 ,fPIDquality(t.fPIDquality)
184 ,fClusterOwner(kTRUE)
185 ,fPIDmethod(t.fPIDmethod)
186 ,fStopped(t.fStopped)
189 ,fNRotate(t.fNRotate)
191 ,fNExpected(t.fNExpected)
193 ,fNExpectedLast(t.fNExpectedLast)
195 ,fChi2Last(t.fChi2Last)
202 for (Int_t i = 0; i < kNplane; i++) {
203 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
204 fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
206 fTimBinPlane[i] = t.fTimBinPlane[i];
207 fTracklets[i] = t.fTracklets[i];
211 fTrackletIndex[i] = t.fTrackletIndex[i];
214 Int_t n = t.GetNumberOfClusters();
215 SetNumberOfClusters(n);
217 for (Int_t i = 0; i < n; i++) {
218 fIndex[i] = t.fIndex[i];
219 fIndexBackup[i] = t.fIndex[i];
220 fdQdl[i] = t.fdQdl[i];
221 if (fClusterOwner && t.fClusters[i]) {
222 fClusters[i] = new AliTRDcluster(*(t.fClusters[i]));
225 fClusters[i] = t.fClusters[i];
229 for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
236 for (Int_t i = 0; i < 3;i++) {
237 fBudget[i] = t.fBudget[i];
240 for(Int_t ispec = 0; ispec<AliPID::kSPECIES; ispec++) fPID[ispec] = t.fPID[ispec];
243 //_____________________________________________________________________________
244 AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
247 ,fdEdx(t.GetPIDsignal())
250 ,fClusterOwner(kFALSE)
265 // Constructor from AliTPCtrack or AliITStrack
268 SetLabel(t.GetLabel());
270 SetMass(t.GetMass());
271 SetNumberOfClusters(0);
273 for (Int_t i = 0; i < kNplane; i++) {
274 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
275 fdEdxPlane[i][j] = 0.0;
277 fTimBinPlane[i] = -1;
281 fTrackletIndex[i] = -1;
283 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
284 fPID[ispec] = 1.0 / AliPID::kSPECIES;
287 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
294 for (Int_t i = 0; i < 3; i++) {
300 //_____________________________________________________________________________
301 AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
304 ,fdEdx(t.GetTRDsignal())
307 ,fClusterOwner(kFALSE)
322 // Constructor from AliESDtrack
325 SetLabel(t.GetLabel());
327 SetMass(t.GetMass());
328 SetNumberOfClusters(t.GetTRDclusters(fIndex));
330 Int_t ncl = t.GetTRDclusters(fIndexBackup);
331 for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) {
336 for (Int_t i = 0; i < kNplane; i++) {
337 for (Int_t j = 0; j < AliTRDCalPID::kNSlicesLQ; j++) {
338 fdEdxPlane[i][j] = t.GetTRDslice(i,j);
340 fTimBinPlane[i] = t.GetTRDTimBin(i);
344 fTrackletIndex[i] = -1;
347 const AliExternalTrackParam *par = &t;
348 if (t.GetStatus() & AliESDtrack::kTRDbackup) {
349 par = t.GetOuterParam();
351 AliError("No backup info!");
358 ,par->GetCovariance());
360 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
365 for (Int_t i = 0; i < 3; i++) {
368 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
369 fPID[ispec] = t.GetTRDpid(ispec);
372 if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
378 t.GetIntegratedTimes(times);
379 SetIntegratedTimes(times);
380 SetIntegratedLength(t.GetIntegratedLength());
384 //____________________________________________________________________________
385 AliTRDtrack::~AliTRDtrack()
398 while ((icluster < kMAXCLUSTERSPERTRACK) && fClusters[icluster]) {
399 delete fClusters[icluster];
400 fClusters[icluster] = 0x0;
407 //____________________________________________________________________________
408 Float_t AliTRDtrack::StatusForTOF()
411 // Defines the status of the TOF extrapolation
414 // Definition of res ????
415 Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0)))
416 * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0);
417 res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
422 //_____________________________________________________________________________
423 Int_t AliTRDtrack::Compare(const TObject *o) const
426 // Compares tracks according to their Y2 or curvature
429 AliTRDtrack *t = (AliTRDtrack *) o;
431 Double_t co = TMath::Abs(t->GetC());
432 Double_t c = TMath::Abs(GetC());
445 //_____________________________________________________________________________
446 void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
449 // Calculates the truncated dE/dx within the "low" and "up" cuts.
452 // Array to sort the dEdx values according to amplitude
453 Float_t sorted[kMAXCLUSTERSPERTRACK];
456 // Require at least 10 clusters for a dedx measurement
461 // Can fdQdl be negative ????
462 for (Int_t i = 0; i < fNdedx; i++) {
463 sorted[i] = TMath::Abs(fdQdl[i]);
465 // Sort the dedx values by amplitude
466 Int_t *index = new Int_t[fNdedx];
467 TMath::Sort(fNdedx, sorted, index, kFALSE);
469 // Sum up the truncated charge between lower and upper bounds
470 Int_t nl = Int_t(low * fNdedx);
471 Int_t nu = Int_t( up * fNdedx);
472 for (Int_t i = nl; i <= nu; i++) {
473 fdEdx += sorted[index[i]];
475 fdEdx /= (nu - nl + 1.0);
481 //_____________________________________________________________________________
482 void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
485 // Set fdEdxPlane and fTimBinPlane and also get the
486 // Time bin for Max. Cluster
489 // Prashant Shukla (shukla@physi.uni-heidelberg.de)
490 // Alexandru Bercuci (A.Bercuci@gsi.de)
493 // Max charge in chamber
494 Double_t maxcharge[kNplane];
495 // Number of clusters attached to track per chamber and slice
496 Int_t nCluster[kNplane][AliTRDCalPID::kNSlicesLQ];
497 // Number of time bins in chamber
498 Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();
499 Int_t plane; // Plane of current cluster
500 Int_t tb; // Time bin of current cluster
501 Int_t slice; // Current slice
502 AliTRDcluster *cluster = 0x0; // Pointer to current cluster
504 // Reset class and local counters/variables
505 for (Int_t iPlane = 0; iPlane < kNplane; iPlane++) {
506 fTimBinPlane[iPlane] = -1;
507 maxcharge[iPlane] = 0.0;
508 for (Int_t iSlice = 0; iSlice < AliTRDCalPID::kNSlicesLQ; iSlice++) {
509 fdEdxPlane[iPlane][iSlice] = 0.0;
510 nCluster[iPlane][iSlice] = 0;
514 // Start looping over clusters attached to this track
515 for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
517 cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
518 if (!cluster) continue;
520 // Read info from current cluster
521 plane = AliTRDgeometry::GetLayer(cluster->GetDetector());
522 if (plane < 0 || plane >= kNplane) {
523 AliError(Form("Wrong plane %d", plane));
527 tb = cluster->GetLocalTimeBin();
528 if ((tb < 0) || (tb >= ntb)) {
529 AliWarning(Form("time bin < 0 or > %d in cluster %d", ntb, iClus));
530 AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
534 slice = tb * AliTRDCalPID::kNSlicesLQ / ntb;
536 fdEdxPlane[plane][slice] += fdQdl[iClus];
537 if (fdQdl[iClus] > maxcharge[plane]) {
538 maxcharge[plane] = fdQdl[iClus];
539 fTimBinPlane[plane] = tb;
542 nCluster[plane][slice]++;
544 } // End of loop over cluster
546 // Normalize fdEdxPlane to number of clusters and set track segments
547 for (Int_t iPlane = 0; iPlane < kNplane; iPlane++) {
548 for (Int_t iSlice = 0; iSlice < AliTRDCalPID::kNSlicesLQ; iSlice++) {
549 if (nCluster[iPlane][iSlice]) {
550 fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
557 //_____________________________________________________________________________
558 void AliTRDtrack::CookdEdxNN(Float_t *dedx)
561 // This function calcuates the input for the neural networks
562 // which are used for the PID.
565 //number of time bins in chamber
566 Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS();
568 Int_t plane; // plane of current cluster
569 Int_t tb; // time bin of current cluster
570 Int_t slice; // curent slice
571 AliTRDcluster *cluster = 0x0; // pointer to current cluster
572 const Int_t kMLPscale = 16000; // scaling of the MLP input to be smaller than 1
574 // Reset class and local contors/variables
575 for (Int_t iPlane = 0; iPlane < kNplane; iPlane++){
576 for (Int_t iSlice = 0; iSlice < AliTRDCalPID::kNSlicesNN; iSlice++) {
577 *(dedx + (AliTRDCalPID::kNSlicesNN * iPlane) + iSlice) = 0.0;
581 // Start looping over clusters attached to this track
582 for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
584 cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
589 // Read info from current cluster
590 plane = AliTRDgeometry::GetLayer(cluster->GetDetector());
591 if (plane < 0 || plane >= kNplane) {
592 AliError(Form("Wrong plane %d",plane));
596 tb = cluster->GetLocalTimeBin();
597 if (tb == 0 || tb >= ntb) {
598 AliWarning(Form("time bin 0 or > %d in cluster %d",ntb,iClus));
599 AliInfo(Form("dQ/dl %f",fdQdl[iClus]));
603 slice = tb * AliTRDCalPID::kNSlicesNN / ntb;
605 *(dedx+(AliTRDCalPID::kNSlicesNN * plane) + slice) += fdQdl[iClus]/kMLPscale;
607 } // End of loop over cluster
611 //_____________________________________________________________________________
612 void AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
615 // Set the momenta for a track segment in a given plane
620 AliError(Form("Trying to access out of range plane (%d)", plane));
624 fSnp[plane] = GetSnp();
625 fTgl[plane] = GetTgl();
628 fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
632 //_____________________________________________________________________________
633 Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
636 // Calculate the track length of a track segment in a given plane
640 (plane >= kNplane)) {
644 return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
645 / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane])
646 / (1.0 + fTgl[plane]*fTgl[plane]));
650 //_____________________________________________________________________________
651 Bool_t AliTRDtrack::CookPID(Int_t &pidQuality)
654 // This function calculates the PID probabilities based on TRD signals
656 // The method produces probabilities based on the charge
657 // and the position of the maximum time bin in each layer.
658 // The dE/dx information can be used as global charge or 2 to 3
659 // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
663 // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
666 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
668 AliError("No access to calibration data");
672 // Retrieve the CDB container class with the probability distributions
673 const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDpidUtil::kLQ);
675 AliError("No access to AliTRDCalPID");
679 // Calculate the input for the NN if fPIDmethod is kNN
680 Float_t ldEdxNN[AliTRDgeometry::kNlayer * AliTRDCalPID::kNSlicesNN], *dedx = 0x0;
681 if(fPIDmethod == kNN) {
682 CookdEdxNN(&ldEdxNN[0]);
685 // Sets the a priori probabilities
686 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
687 fPID[ispec] = 1.0 / AliPID::kSPECIES;
690 if(AliPID::kSPECIES != 5){
691 AliError("Probabilities array defined only for 5 values. Please modify !!");
696 Float_t length; // track segment length in chamber
698 // Skip tracks which have no TRD signal at all
699 if (fdEdx == 0.) return kFALSE;
701 for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNlayer; iPlane++) {
703 length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
704 / TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane])
705 / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
708 if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0.
709 || fTimBinPlane[iPlane] == -1.) continue;
711 // this track segment has fulfilled all requierments
714 // Get the probabilities for the different particle species
715 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
718 dedx = fdEdxPlane[iPlane];
721 dedx = &ldEdxNN[iPlane*AliTRDCalPID::kNSlicesNN];
724 fPID[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane);
729 if (pidQuality == 0) {
733 // normalize probabilities
734 Double_t probTotal = 0.0;
735 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
736 probTotal += fPID[iSpecies];
739 if (probTotal <= 0.0) {
740 AliWarning("The total probability over all species <= 0.");
741 AliWarning("This may be caused by some error in the reference histograms.");
745 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
746 fPID[iSpecies] /= probTotal;
753 //_____________________________________________________________________________
754 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
757 // Propagates this track to a reference plane defined by "xk" [cm]
758 // correcting for the mean crossed material.
760 // "xx0" - thickness/rad.length [units of the radiation length]
761 // "xrho" - thickness*density [g/cm^2]
768 Double_t oldX = GetX();
769 Double_t oldY = GetY();
770 Double_t oldZ = GetZ();
772 Double_t bz = GetBz();
774 if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
784 if (IsStartedTimeIntegral()) {
785 Double_t l2 = TMath::Sqrt((x-oldX)*(x-oldX)
787 + (z-oldZ)*(z-oldZ));
788 Double_t crv = GetC();
789 if (TMath::Abs(l2*crv) > 0.0001) {
790 // Make correction for curvature if neccesary
791 l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX)
792 + (y-oldY)*(y-oldY));
793 l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
794 l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
800 if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
807 Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
808 Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
809 if ((beta2 < 1.0e-10) ||
810 ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
814 Double_t dE = 0.153e-3 / beta2
815 * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
820 // Suspicious part - think about it ?
821 Double_t kinE = TMath::Sqrt(p2);
822 if (dE > 0.8*kinE) dE = 0.8 * kinE; //
823 if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch
830 Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation
831 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
833 fCee += fX*fX * sigmac2;
842 //_____________________________________________________________________________
843 Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq
844 , Int_t index, Double_t h01)
847 // Assignes the found cluster <c> to the track and updates
848 // track information.
849 // <chisq> : predicted chi2
851 // <h01> : Tilting factor
854 Double_t p[2] = { c->GetY()
856 Double_t sy2 = c->GetSigmaY2() * 4.0;
857 Double_t sz2 = c->GetSigmaZ2() * 4.0;
858 Double_t cov[3] = { sy2 + h01*h01*sz2
860 , sz2 + h01*h01*sy2 };
862 if (!AliExternalTrackParam::Update(p,cov)) {
866 Int_t n = GetNumberOfClusters();
868 SetNumberOfClusters(n+1);
870 SetChi2(GetChi2()+chisq);
876 //_____________________________________________________________________________
877 Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
878 , Double_t h01, Int_t /*plane*/, Int_t /*tid*/)
881 // Assignes the found cluster <c> to the track and
882 // updates track information
883 // <chisq> : predicted chi2
885 // <h01> : Tilting factor
887 // Difference to Update(AliTRDcluster *c): cluster is added to fClusters
890 Double_t p[2] = { c->GetY()
892 Double_t sy2 = c->GetSigmaY2() * 4.0;
893 Double_t sz2 = c->GetSigmaZ2() * 4.0;
894 Double_t cov[3] = { sy2 + h01*h01*sz2
896 , sz2 + h01*h01*sy2 };
898 if (!AliExternalTrackParam::Update(p,cov)) {
902 AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
904 // Register cluster to track
905 Int_t n = GetNumberOfClusters();
908 SetNumberOfClusters(n+1);
909 SetChi2(GetChi2() + chisq);
915 //_____________________________________________________________________________
916 Bool_t AliTRDtrack::Update(const AliTRDtracklet &t, Double_t chisq, Int_t index)
919 // Assignes the found tracklet <t> to the track
920 // and updates track information
921 // <chisq> : predicted chi2
925 Double_t h01 = t.GetTilt(); // Is this correct????
926 Double_t p[2] = { t.GetY(), t.GetZ() };
927 Double_t sy2 = t.GetTrackletSigma2(); // Error pad-column
928 Double_t sz2 = 100000.0; // Error pad-row (????)
929 Double_t cov[3] = { sy2 + h01*h01*sz2 // Does this have any sense????
931 , sz2 + h01*h01*sy2 };
932 if (!AliExternalTrackParam::Update(p,cov)) {
936 Int_t n = GetNumberOfClusters();
938 SetNumberOfClusters(n+1);
939 SetChi2(GetChi2()+chisq);
945 //_____________________________________________________________________________
946 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
949 // Rotates track parameters in R*phi plane
950 // if absolute rotation alpha is in global system
951 // otherwise alpha rotation is relative to the current rotation angle
961 return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
965 //_____________________________________________________________________________
966 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
969 // Returns the track chi2
972 Double_t p[2] = { c->GetY()
974 Double_t sy2 = c->GetSigmaY2() * 4.0;
975 Double_t sz2 = c->GetSigmaZ2() * 4.0;
976 Double_t cov[3] = { sy2 + h01*h01*sz2
978 , sz2 + h01*h01*sy2 };
980 return AliExternalTrackParam::GetPredictedChi2(p,cov);
984 //_____________________________________________________________________________
985 void AliTRDtrack::MakeBackupTrack()
988 // Creates a backup track
994 fBackupTrack = new AliTRDtrack(*this);
998 //_____________________________________________________________________________
999 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
1002 // Find a prolongation at given x
1003 // Return 0 if it does not exist
1006 Double_t bz = GetBz();
1008 if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
1011 if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
1019 //_____________________________________________________________________________
1020 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
1023 // Propagate track to given x position
1024 // Works inside of the 20 degree segmentation
1025 // (local cooordinate frame for TRD , TPC, TOF)
1027 // Material budget from geo manager
1035 const Double_t kAlphac = TMath::Pi()/9.0;
1036 const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
1038 // Critical alpha - cross sector indication
1039 Double_t dir = (GetX() > xr) ? -1.0 : 1.0;
1042 for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1045 GetProlongation(x,y,z);
1046 xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1047 xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
1050 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1052 if ((param[0] > 0) &&
1054 PropagateTo(x,param[1],param[0]*param[4]);
1057 if (GetY() > GetX()*kTalphac) {
1060 if (GetY() < -GetX()*kTalphac) {
1072 //_____________________________________________________________________________
1073 Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step)
1076 // Propagate track to the radial position
1077 // Rotation always connected to the last track position
1085 Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
1087 Double_t dir = (radius > r) ? -1.0 : 1.0;
1089 for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1092 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1093 Rotate(alpha,kTRUE);
1095 GetProlongation(x,y,z);
1096 xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1097 xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1100 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1101 if (param[1] <= 0) {
1102 param[1] = 100000000;
1104 PropagateTo(x,param[1],param[0]*param[4]);
1109 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1110 Rotate(alpha,kTRUE);
1112 GetProlongation(r,y,z);
1113 xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1114 xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1117 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1119 if (param[1] <= 0) {
1120 param[1] = 100000000;
1122 PropagateTo(r,param[1],param[0]*param[4]);
1128 //_____________________________________________________________________________
1129 Int_t AliTRDtrack::GetSector() const
1132 // Return the current sector
1135 return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
1136 % AliTRDgeometry::kNsector;
1140 //_____________________________________________________________________________
1141 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
1144 // The sampled energy loss
1147 Double_t s = GetSnp();
1148 Double_t t = GetTgl();
1149 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1154 //_____________________________________________________________________________
1155 void AliTRDtrack::SetSampledEdx(Float_t q)
1158 // The sampled energy loss
1161 Double_t s = GetSnp();
1162 Double_t t = GetTgl();
1163 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1169 //_____________________________________________________________________________
1170 Double_t AliTRDtrack::GetBz() const
1173 // Returns Bz component of the magnetic field (kG)
1176 if (AliTracker::UniformField()) {
1177 return AliTracker::GetBz();
1182 return AliTracker::GetBz(r);