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"
21 #include "AliESDtrack.h"
23 #include "AliTRDgeometry.h"
24 #include "AliTRDcluster.h"
25 #include "AliTRDtrack.h"
26 #include "AliTRDcalibDB.h"
27 #include "Cal/AliTRDCalPID.h"
31 ///////////////////////////////////////////////////////////////////////////////
33 // Represents a reconstructed TRD track //
34 // Local TRD Kalman track //
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 < kNslice; j++) {
65 fdEdxPlane[i][j] = 0.0;
73 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
80 for (Int_t i = 0; i < 3; i++) {
86 //_____________________________________________________________________________
87 AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index
88 , const Double_t p[5], const Double_t cov[15]
89 , Double_t x, Double_t alpha)
95 ,fClusterOwner(kFALSE)
110 // The main AliTRDtrack constructor.
113 Double_t cnv = 1.0 / (GetBz() * kB2C);
115 Double_t pp[5] = { p[0]
121 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
122 Double_t c32 = x*cov[13] - cov[ 8];
123 Double_t c20 = x*cov[10] - cov[ 3];
124 Double_t c21 = x*cov[11] - cov[ 4];
125 Double_t c42 = x*cov[14] - cov[12];
127 Double_t cc[15] = { cov[ 0]
130 , cov[ 6], cov[ 7], c32, cov[ 9]
131 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
134 SetNumberOfClusters(1);
138 for (Int_t i = 0; i < kNplane; i++) {
139 for (Int_t j = 0; j < kNslice; j++) {
140 fdEdxPlane[i][j] = 0.0;
142 fTimBinPlane[i] = -1;
148 Double_t q = TMath::Abs(c->GetQ());
149 Double_t s = GetSnp();
150 Double_t t = GetTgl();
152 q *= TMath::Sqrt((1-s*s)/(1+t*t));
156 for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
163 for (Int_t i = 0; i < 3;i++) {
169 //_____________________________________________________________________________
170 AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
172 ,fSeedLab(t.GetSeedLabel())
175 ,fPIDquality(t.fPIDquality)
176 ,fClusterOwner(kTRUE)
177 ,fPIDmethod(t.fPIDmethod)
178 ,fStopped(t.fStopped)
181 ,fNRotate(t.fNRotate)
183 ,fNExpected(t.fNExpected)
185 ,fNExpectedLast(t.fNExpectedLast)
187 ,fChi2Last(t.fChi2Last)
194 for (Int_t i = 0; i < kNplane; i++) {
195 for (Int_t j = 0; j < kNslice; j++) {
196 fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
198 fTimBinPlane[i] = t.fTimBinPlane[i];
199 fTracklets[i] = t.fTracklets[i];
205 Int_t n = t.GetNumberOfClusters();
206 SetNumberOfClusters(n);
208 for (Int_t i = 0; i < n; i++) {
209 fIndex[i] = t.fIndex[i];
210 fIndexBackup[i] = t.fIndex[i];
211 fdQdl[i] = t.fdQdl[i];
212 if (fClusterOwner && t.fClusters[i]) {
213 fClusters[i] = new AliTRDcluster(*(t.fClusters[i]));
216 fClusters[i] = t.fClusters[i];
220 for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
227 for (Int_t i = 0; i < 3;i++) {
228 fBudget[i] = t.fBudget[i];
231 for(Int_t ispec = 0; ispec<AliPID::kSPECIES; ispec++) fPID[ispec] = t.fPID[ispec];
234 //_____________________________________________________________________________
235 AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
238 ,fdEdx(t.GetPIDsignal())
241 ,fClusterOwner(kFALSE)
256 // Constructor from AliTPCtrack or AliITStrack
259 SetLabel(t.GetLabel());
261 SetMass(t.GetMass());
262 SetNumberOfClusters(0);
264 for (Int_t i = 0; i < kNplane; i++) {
265 for (Int_t j = 0; j < kNslice; j++) {
266 fdEdxPlane[i][j] = 0.0;
268 fTimBinPlane[i] = -1;
274 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
281 for (Int_t i = 0; i < 3; i++) {
287 //_____________________________________________________________________________
288 AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
291 ,fdEdx(t.GetTRDsignal())
294 ,fClusterOwner(kFALSE)
309 // Constructor from AliESDtrack
312 SetLabel(t.GetLabel());
314 SetMass(t.GetMass());
315 SetNumberOfClusters(t.GetTRDclusters(fIndex));
317 Int_t ncl = t.GetTRDclusters(fIndexBackup);
318 for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) {
323 for (Int_t i = 0; i < kNplane; i++) {
324 for (Int_t j = 0; j < kNslice; j++) {
325 fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
327 fTimBinPlane[i] = t.GetTRDTimBin(i);
333 const AliExternalTrackParam *par = &t;
334 if (t.GetStatus() & AliESDtrack::kTRDbackup) {
335 par = t.GetOuterParam();
337 AliError("No backup info!");
344 ,par->GetCovariance());
346 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
351 for (Int_t i = 0; i < 3; i++) {
355 if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
361 t.GetIntegratedTimes(times);
362 SetIntegratedTimes(times);
363 SetIntegratedLength(t.GetIntegratedLength());
367 //____________________________________________________________________________
368 AliTRDtrack::~AliTRDtrack()
381 while ((icluster < kMAXCLUSTERSPERTRACK) && fClusters[icluster]) {
382 delete fClusters[icluster];
383 fClusters[icluster] = 0x0;
390 //____________________________________________________________________________
391 Float_t AliTRDtrack::StatusForTOF()
394 // Defines the status of the TOF extrapolation
397 // Definition of res ????
398 Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0)))
399 * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0);
400 res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
405 //_____________________________________________________________________________
406 Int_t AliTRDtrack::Compare(const TObject *o) const
409 // Compares tracks according to their Y2 or curvature
412 AliTRDtrack *t = (AliTRDtrack *) o;
414 Double_t co = TMath::Abs(t->GetC());
415 Double_t c = TMath::Abs(GetC());
428 //_____________________________________________________________________________
429 void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
432 // Calculates the truncated dE/dx within the "low" and "up" cuts.
435 // Array to sort the dEdx values according to amplitude
436 Float_t sorted[kMAXCLUSTERSPERTRACK];
439 // Require at least 10 clusters for a dedx measurement
444 // Can fdQdl be negative ????
445 for (Int_t i = 0; i < fNdedx; i++) {
446 sorted[i] = TMath::Abs(fdQdl[i]);
448 // Sort the dedx values by amplitude
449 Int_t *index = new Int_t[fNdedx];
450 TMath::Sort(fNdedx, sorted, index, kFALSE);
452 // Sum up the truncated charge between lower and upper bounds
453 Int_t nl = Int_t(low * fNdedx);
454 Int_t nu = Int_t( up * fNdedx);
455 for (Int_t i = nl; i <= nu; i++) {
456 fdEdx += sorted[index[i]];
458 fdEdx /= (nu - nl + 1.0);
464 //_____________________________________________________________________________
465 void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
468 // Set fdEdxPlane and fTimBinPlane and also get the
469 // Time bin for Max. Cluster
472 // Prashant Shukla (shukla@physi.uni-heidelberg.de)
473 // Alexandru Bercuci (A.Bercuci@gsi.de)
476 // Max charge in chamber
477 Double_t maxcharge[AliESDtrack::kNPlane];
478 // Number of clusters attached to track per chamber and slice
479 Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
480 // Number of time bins in chamber
481 Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
482 Int_t plane; // Plane of current cluster
483 Int_t tb; // Time bin of current cluster
484 Int_t slice; // Current slice
485 AliTRDcluster *cluster = 0x0; // Pointer to current cluster
487 // Reset class and local counters/variables
488 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
489 fTimBinPlane[iPlane] = -1;
490 maxcharge[iPlane] = 0.0;
491 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
492 fdEdxPlane[iPlane][iSlice] = 0.0;
493 nCluster[iPlane][iSlice] = 0;
497 // Start looping over clusters attached to this track
498 for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
500 cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
501 if (!cluster) continue;
503 // Read info from current cluster
504 plane = AliTRDgeometry::GetPlane(cluster->GetDetector());
505 if (plane < 0 || plane >= AliESDtrack::kNPlane) {
506 AliError(Form("Wrong plane %d", plane));
510 tb = cluster->GetLocalTimeBin();
511 if ((tb < 0) || (tb >= ntb)) {
512 AliWarning(Form("time bin < 0 or > %d in cluster %d", ntb, iClus));
513 AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
517 slice = tb * AliESDtrack::kNSlice / ntb;
519 fdEdxPlane[plane][slice] += fdQdl[iClus];
520 if (fdQdl[iClus] > maxcharge[plane]) {
521 maxcharge[plane] = fdQdl[iClus];
522 fTimBinPlane[plane] = tb;
525 nCluster[plane][slice]++;
527 } // End of loop over cluster
529 // Normalize fdEdxPlane to number of clusters and set track segments
530 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
531 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
532 if (nCluster[iPlane][iSlice]) {
533 fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
540 //_____________________________________________________________________________
541 void AliTRDtrack::CookdEdxNN(Float_t *dedx)
544 // This function calcuates the input for the neural networks
545 // which are used for the PID.
548 //number of time bins in chamber
549 Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
551 Int_t plane; // plane of current cluster
552 Int_t tb; // time bin of current cluster
553 Int_t slice; // curent slice
554 AliTRDcluster *cluster = 0x0; // pointer to current cluster
555 const Int_t kMLPscale = 16000; // scaling of the MLP input to be smaller than 1
557 // Reset class and local contors/variables
558 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++){
559 for (Int_t iSlice = 0; iSlice < kNMLPslice; iSlice++) {
560 *(dedx + (kNMLPslice * iPlane) + iSlice) = 0.0;
564 // Start looping over clusters attached to this track
565 for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
567 cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
572 // Read info from current cluster
573 plane = AliTRDgeometry::GetPlane(cluster->GetDetector());
574 if (plane < 0 || plane >= AliESDtrack::kNPlane) {
575 AliError(Form("Wrong plane %d",plane));
579 tb = cluster->GetLocalTimeBin();
580 if (tb == 0 || tb >= ntb) {
581 AliWarning(Form("time bin 0 or > %d in cluster %d",ntb,iClus));
582 AliInfo(Form("dQ/dl %f",fdQdl[iClus]));
586 slice = tb * kNMLPslice / ntb;
588 *(dedx+(kNMLPslice * plane) + slice) += fdQdl[iClus]/kMLPscale;
590 } // End of loop over cluster
594 //_____________________________________________________________________________
595 void AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
598 // Set the momenta for a track segment in a given plane
603 AliError(Form("Trying to access out of range plane (%d)", plane));
607 fSnp[plane] = GetSnp();
608 fTgl[plane] = GetTgl();
611 fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
615 //_____________________________________________________________________________
616 Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
619 // Calculate the track length of a track segment in a given plane
623 (plane >= kNplane)) {
627 return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
628 / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane])
629 / (1.0 + fTgl[plane]*fTgl[plane]));
633 //_____________________________________________________________________________
634 Bool_t AliTRDtrack::CookPID(Int_t &pidQuality)
637 // This function calculates the PID probabilities based on TRD signals
639 // The method produces probabilities based on the charge
640 // and the position of the maximum time bin in each layer.
641 // The dE/dx information can be used as global charge or 2 to 3
642 // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
646 // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
649 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
651 AliError("No access to calibration data");
655 // Retrieve the CDB container class with the probability distributions
656 const AliTRDCalPID *pd = calibration->GetPIDObject(fPIDmethod == kNN ? 0 : 1);
658 AliError("No access to AliTRDCalPID");
662 // Calculate the input for the NN if fPIDmethod is kNN
663 Float_t ldEdxNN[AliTRDCalPID::kNPlane * kNMLPslice], *dedx = 0x0;
664 if(fPIDmethod == kNN) {
665 CookdEdxNN(&ldEdxNN[0]);
668 // Sets the a priori probabilities
669 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
670 fPID[ispec] = 1.0 / AliPID::kSPECIES;
673 if(AliPID::kSPECIES != 5){
674 AliError("Probabilities array defined only for 5 values. Please modify !!");
679 Float_t length; // track segment length in chamber
681 // Skip tracks which have no TRD signal at all
682 if (fdEdx == 0.) return kFALSE;
684 for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
686 length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
687 / TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane])
688 / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
691 if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0.
692 || fTimBinPlane[iPlane] == -1.) continue;
694 // this track segment has fulfilled all requierments
697 // Get the probabilities for the different particle species
698 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
701 dedx = fdEdxPlane[iPlane];
704 dedx = &ldEdxNN[iPlane*kNMLPslice];
707 fPID[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane);
712 if (pidQuality == 0) {
716 // normalize probabilities
717 Double_t probTotal = 0.0;
718 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
719 probTotal += fPID[iSpecies];
722 if (probTotal <= 0.0) {
723 AliWarning("The total probability over all species <= 0.");
724 AliWarning("This may be caused by some error in the reference histograms.");
728 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
729 fPID[iSpecies] /= probTotal;
736 //_____________________________________________________________________________
737 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
740 // Propagates this track to a reference plane defined by "xk" [cm]
741 // correcting for the mean crossed material.
743 // "xx0" - thickness/rad.length [units of the radiation length]
744 // "xrho" - thickness*density [g/cm^2]
751 Double_t oldX = GetX();
752 Double_t oldY = GetY();
753 Double_t oldZ = GetZ();
755 Double_t bz = GetBz();
757 if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
767 if (IsStartedTimeIntegral()) {
768 Double_t l2 = TMath::Sqrt((x-oldX)*(x-oldX)
770 + (z-oldZ)*(z-oldZ));
771 Double_t crv = GetC();
772 if (TMath::Abs(l2*crv) > 0.0001) {
773 // Make correction for curvature if neccesary
774 l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX)
775 + (y-oldY)*(y-oldY));
776 l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
777 l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
783 if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
790 Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
791 Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
792 if ((beta2 < 1.0e-10) ||
793 ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
797 Double_t dE = 0.153e-3 / beta2
798 * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
803 // Suspicious part - think about it ?
804 Double_t kinE = TMath::Sqrt(p2);
805 if (dE > 0.8*kinE) dE = 0.8 * kinE; //
806 if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch
813 Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation
814 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
816 fCee += fX*fX * sigmac2;
825 //_____________________________________________________________________________
826 Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq
827 , Int_t index, Double_t h01)
830 // Assignes the found cluster <c> to the track and updates
831 // track information.
832 // <chisq> : predicted chi2
834 // <h01> : Tilting factor
837 Double_t p[2] = { c->GetY()
839 Double_t sy2 = c->GetSigmaY2() * 4.0;
840 Double_t sz2 = c->GetSigmaZ2() * 4.0;
841 Double_t cov[3] = { sy2 + h01*h01*sz2
843 , sz2 + h01*h01*sy2 };
845 if (!AliExternalTrackParam::Update(p,cov)) {
849 Int_t n = GetNumberOfClusters();
851 SetNumberOfClusters(n+1);
853 SetChi2(GetChi2()+chisq);
859 //_____________________________________________________________________________
860 Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
861 , Double_t h01, Int_t /*plane*/, Int_t /*tid*/)
864 // Assignes the found cluster <c> to the track and
865 // updates track information
866 // <chisq> : predicted chi2
868 // <h01> : Tilting factor
870 // Difference to Update(AliTRDcluster *c): cluster is added to fClusters
873 Double_t p[2] = { c->GetY()
875 Double_t sy2 = c->GetSigmaY2() * 4.0;
876 Double_t sz2 = c->GetSigmaZ2() * 4.0;
877 Double_t cov[3] = { sy2 + h01*h01*sz2
879 , sz2 + h01*h01*sy2 };
881 if (!AliExternalTrackParam::Update(p,cov)) {
885 AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
887 // Register cluster to track
888 Int_t n = GetNumberOfClusters();
891 SetNumberOfClusters(n+1);
892 SetChi2(GetChi2() + chisq);
898 //_____________________________________________________________________________
899 Bool_t AliTRDtrack::Update(const AliTRDtracklet &t, Double_t chisq, Int_t index)
902 // Assignes the found tracklet <t> to the track
903 // and updates track information
904 // <chisq> : predicted chi2
908 Double_t h01 = t.GetTilt(); // Is this correct????
909 Double_t p[2] = { t.GetY(), t.GetZ() };
910 Double_t sy2 = t.GetTrackletSigma2(); // Error pad-column
911 Double_t sz2 = 100000.0; // Error pad-row (????)
912 Double_t cov[3] = { sy2 + h01*h01*sz2 // Does this have any sense????
914 , sz2 + h01*h01*sy2 };
915 if (!AliExternalTrackParam::Update(p,cov)) {
919 Int_t n = GetNumberOfClusters();
921 SetNumberOfClusters(n+1);
922 SetChi2(GetChi2()+chisq);
928 //_____________________________________________________________________________
929 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
932 // Rotates track parameters in R*phi plane
933 // if absolute rotation alpha is in global system
934 // otherwise alpha rotation is relative to the current rotation angle
944 return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
948 //_____________________________________________________________________________
949 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
952 // Returns the track chi2
955 Double_t p[2] = { c->GetY()
957 Double_t sy2 = c->GetSigmaY2() * 4.0;
958 Double_t sz2 = c->GetSigmaZ2() * 4.0;
959 Double_t cov[3] = { sy2 + h01*h01*sz2
961 , sz2 + h01*h01*sy2 };
963 return AliExternalTrackParam::GetPredictedChi2(p,cov);
967 //_____________________________________________________________________________
968 void AliTRDtrack::MakeBackupTrack()
971 // Creates a backup track
977 fBackupTrack = new AliTRDtrack(*this);
981 //_____________________________________________________________________________
982 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
985 // Find a prolongation at given x
986 // Return 0 if it does not exist
989 Double_t bz = GetBz();
991 if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
994 if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
1002 //_____________________________________________________________________________
1003 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
1006 // Propagate track to given x position
1007 // Works inside of the 20 degree segmentation
1008 // (local cooordinate frame for TRD , TPC, TOF)
1010 // Material budget from geo manager
1018 const Double_t kAlphac = TMath::Pi()/9.0;
1019 const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
1021 // Critical alpha - cross sector indication
1022 Double_t dir = (GetX() > xr) ? -1.0 : 1.0;
1025 for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1028 GetProlongation(x,y,z);
1029 xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1030 xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
1033 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1035 if ((param[0] > 0) &&
1037 PropagateTo(x,param[1],param[0]*param[4]);
1040 if (GetY() > GetX()*kTalphac) {
1043 if (GetY() < -GetX()*kTalphac) {
1055 //_____________________________________________________________________________
1056 Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step)
1059 // Propagate track to the radial position
1060 // Rotation always connected to the last track position
1068 Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
1070 Double_t dir = (radius > r) ? -1.0 : 1.0;
1072 for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1075 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1076 Rotate(alpha,kTRUE);
1078 GetProlongation(x,y,z);
1079 xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1080 xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1083 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1084 if (param[1] <= 0) {
1085 param[1] = 100000000;
1087 PropagateTo(x,param[1],param[0]*param[4]);
1092 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1093 Rotate(alpha,kTRUE);
1095 GetProlongation(r,y,z);
1096 xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1097 xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1100 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1102 if (param[1] <= 0) {
1103 param[1] = 100000000;
1105 PropagateTo(r,param[1],param[0]*param[4]);
1111 //_____________________________________________________________________________
1112 Int_t AliTRDtrack::GetSector() const
1115 // Return the current sector
1118 return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
1119 % AliTRDgeometry::kNsect;
1123 //_____________________________________________________________________________
1124 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
1127 // The sampled energy loss
1130 Double_t s = GetSnp();
1131 Double_t t = GetTgl();
1132 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1137 //_____________________________________________________________________________
1138 void AliTRDtrack::SetSampledEdx(Float_t q)
1141 // The sampled energy loss
1144 Double_t s = GetSnp();
1145 Double_t t = GetTgl();
1146 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1152 //_____________________________________________________________________________
1153 Double_t AliTRDtrack::GetBz() const
1156 // Returns Bz component of the magnetic field (kG)
1159 if (AliTracker::UniformField()) {
1160 return AliTracker::GetBz();
1165 return AliTracker::GetBz(r);