LCOV - code coverage report
Current view: top level - TRD/TRDbase - AliTRDtrackV1.cxx (source / functions) Hit Total Coverage
Test: coverage.info Lines: 267 493 54.2 %
Date: 2016-06-14 17:26:59 Functions: 27 39 69.2 %

          Line data    Source code
       1             : 
       2             : /**************************************************************************
       3             :  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
       4             :  *                                                                        *
       5             :  * Author: The ALICE Off-line Project.                                    *
       6             :  * Contributors are mentioned in the code where appropriate.              *
       7             :  *                                                                        *
       8             :  * Permission to use, copy, modify and distribute this software and its   *
       9             :  * documentation strictly for non-commercial purposes is hereby granted   *
      10             :  * without fee, provided that the above copyright notice appears in all   *
      11             :  * copies and that both the copyright notice and this permission notice   *
      12             :  * appear in the supporting documentation. The authors make no claims     *
      13             :  * about the suitability of this software for any purpose. It is          *
      14             :  * provided "as is" without express or implied warranty.                  *
      15             :  **************************************************************************/
      16             : 
      17             : /* $Id$ */
      18             : 
      19             : #include "TVectorT.h"
      20             : #include "AliLog.h"
      21             : #include "AliESDtrack.h"
      22             : #include "AliTracker.h"
      23             : 
      24             : #include "AliTRDtrackV1.h"
      25             : #include "AliTRDcluster.h"
      26             : #include "AliTRDcalibDB.h"
      27             : #include "AliTRDReconstructor.h"
      28             : #include "AliTRDPIDResponse.h"
      29             : #include "AliTRDrecoParam.h"
      30             : #include "AliTRDdEdxBaseUtils.h"
      31             : #include "AliTRDdEdxCalibHistArray.h"
      32             : #include "AliTRDdEdxCalibUtils.h"
      33             : #include "AliTRDdEdxReconUtils.h"
      34             : 
      35          48 : ClassImp(AliTRDtrackV1)
      36             : 
      37             : ///////////////////////////////////////////////////////////////////////////////
      38             : //                                                                           //
      39             : //  Represents a reconstructed TRD track                                     //
      40             : //  Local TRD Kalman track                                                   //
      41             : //                                                                           //
      42             : //  Authors:                                                                 //
      43             : //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
      44             : //    Markus Fasel <M.Fasel@gsi.de>                                          //
      45             : //                                                                           //
      46             : ///////////////////////////////////////////////////////////////////////////////
      47             : 
      48             : //_______________________________________________________________
      49          16 : AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
      50          16 :   ,fStatus(0)
      51          16 :   ,fESDid(0)
      52          16 :   ,fDE(0.)
      53          16 :   ,fTruncatedMean(0)
      54          16 :   ,fNchamberdEdx(0)
      55          16 :   ,fNclusterdEdx(0)
      56          16 :   ,fNdEdxSlices(0)
      57          16 :   ,fkReconstructor(NULL)
      58          16 :   ,fBackupTrack(NULL)
      59          16 :   ,fTrackLow(NULL)
      60          16 :   ,fTrackHigh(NULL)
      61          80 : {
      62             :   //
      63             :   // Default constructor
      64             :   //
      65             :   //printf("AliTRDtrackV1::AliTRDtrackV1()\n");
      66             : 
      67         128 :   for(int i =0; i<3; i++) fBudget[i] = 0.;
      68             :   
      69             :   Float_t pid = 1./AliPID::kSPECIES;
      70         192 :   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
      71             : 
      72         224 :   for(int ip=0; ip<kNplane; ip++){
      73          96 :     fTrackletIndex[ip] = -1;
      74          96 :     fTracklet[ip]      = NULL;
      75             :   }
      76          16 :   SetLabel(-123456789); // reset label
      77          32 : }
      78             : 
      79             : //_______________________________________________________________
      80         104 : AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
      81         104 :   ,fStatus(ref.fStatus)
      82         104 :   ,fESDid(ref.fESDid)
      83         104 :   ,fDE(ref.fDE)
      84         104 :   ,fTruncatedMean(ref.fTruncatedMean)
      85         104 :   ,fNchamberdEdx(ref.fNchamberdEdx)
      86         104 :   ,fNclusterdEdx(ref.fNclusterdEdx)
      87         104 :   ,fNdEdxSlices(ref.fNdEdxSlices)
      88         104 :   ,fkReconstructor(ref.fkReconstructor)
      89         104 :   ,fBackupTrack(NULL)
      90         104 :   ,fTrackLow(NULL)
      91         104 :   ,fTrackHigh(NULL)
      92         520 : {
      93             :   //
      94             :   // Copy constructor
      95             :   //
      96             : 
      97             :   //printf("AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &)\n");
      98         104 :   SetBit(kOwner, kFALSE);
      99        1456 :   for(int ip=0; ip<kNplane; ip++){ 
     100         624 :     fTrackletIndex[ip] = ref.fTrackletIndex[ip];
     101         624 :     fTracklet[ip]      = ref.fTracklet[ip];
     102             :   }
     103         230 :   if(ref.fTrackLow) fTrackLow = new AliExternalTrackParam(*ref.fTrackLow);
     104         104 :   if(ref.fTrackHigh) fTrackHigh = new AliExternalTrackParam(*ref.fTrackHigh);
     105             :  
     106         832 :   for (Int_t i = 0; i < 3;i++) fBudget[i]      = ref.fBudget[i];
     107             : 
     108        1248 :         for(Int_t is = 0; is<AliPID::kSPECIES; is++) fPID[is] = ref.fPID[is];
     109             :   
     110         208 :   AliKalmanTrack::SetNumberOfClusters(ref.GetNumberOfClusters());
     111         208 : }
     112             : 
     113             : //_______________________________________________________________
     114         286 : AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
     115         286 :   ,fStatus(0)
     116         286 :   ,fESDid(0)
     117         286 :   ,fDE(0.)
     118         286 :   ,fTruncatedMean(0)
     119         286 :   ,fNchamberdEdx(0)                                                 
     120         286 :   ,fNclusterdEdx(0)
     121         286 :   ,fNdEdxSlices(0)
     122         286 :   ,fkReconstructor(NULL)
     123         286 :   ,fBackupTrack(NULL)
     124         286 :   ,fTrackLow(NULL)
     125         286 :   ,fTrackHigh(NULL)
     126        1430 : {
     127             :   //
     128             :   // Constructor from AliESDtrack
     129             :   //
     130             : 
     131         572 :   SetESDid(t.GetID());
     132         572 :   SetLabel(t.GetLabel());
     133         286 :   SetChi2(0.0);
     134             : 
     135         572 :   SetMass(t.GetMassForTracking());
     136         572 :   AliKalmanTrack::SetNumberOfClusters(t.GetTRDncls()); 
     137         286 :   Int_t ti[]={-1, -1, -1, -1, -1, -1}; t.GetTRDtracklets(&ti[0]);
     138        4004 :   for(int ip=0; ip<kNplane; ip++){ 
     139        1716 :     fTrackletIndex[ip] = ti[ip];
     140        1716 :     fTracklet[ip]      = NULL;
     141             :   }
     142        2288 :   for(int i =0; i<3; i++) fBudget[i] = 0.;
     143             :   
     144             :   Float_t pid = 1./AliPID::kSPECIES;
     145        3432 :   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
     146             : 
     147         286 :   const AliExternalTrackParam *par = &t;
     148             :   /* RS
     149             :   if (t.GetStatus() & AliESDtrack::kTRDbackup) {
     150             :     par = t.GetOuterParam();
     151             :     if (!par) {
     152             :       AliError("No backup info!"); 
     153             :       par = &t;
     154             :     }
     155             :   }
     156             :   */
     157         572 :   Set(par->GetX() 
     158         286 :      ,par->GetAlpha()
     159         286 :      ,par->GetParameter()
     160         286 :      ,par->GetCovariance());
     161             : 
     162         572 :   if(t.GetStatus() & AliESDtrack::kTIME) {
     163         286 :     StartTimeIntegral();
     164         286 :     Double_t times[AliPID::kSPECIESC]; 
     165         286 :     t.GetIntegratedTimes(times,AliPID::kSPECIESC); 
     166         286 :     SetIntegratedTimes(times);
     167         572 :     SetIntegratedLength(t.GetIntegratedLength());
     168         286 :   }
     169         572 : }
     170             : 
     171             : //_______________________________________________________________
     172             : AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 * const trklts, const Double_t p[5], const Double_t cov[15]
     173           0 :              , Double_t x, Double_t alpha) : AliKalmanTrack()
     174           0 :   ,fStatus(0)
     175           0 :   ,fESDid(0)
     176           0 :   ,fDE(0.)
     177           0 :   ,fTruncatedMean(0)
     178           0 :   ,fNchamberdEdx(0)                                                 
     179           0 :   ,fNclusterdEdx(0)
     180           0 :   ,fNdEdxSlices(0)
     181           0 :   ,fkReconstructor(NULL)
     182           0 :   ,fBackupTrack(NULL)
     183           0 :   ,fTrackLow(NULL)
     184           0 :   ,fTrackHigh(NULL)
     185           0 : {
     186             :   //
     187             :   // The stand alone tracking constructor
     188             :   // TEMPORARY !!!!!!!!!!!
     189             :   // to check :
     190             :   // 1. covariance matrix
     191             :   // 2. dQdl calculation
     192             :   //
     193             : 
     194           0 :   Double_t b(GetBz());
     195           0 :   Double_t cnv = (TMath::Abs(b) < 1.e-5) ? 1.e5 : 1./GetBz()/kB2C;
     196             :   
     197           0 :   Double_t pp[5] = { p[0]    
     198           0 :                     , p[1]
     199           0 :                     , p[2]
     200           0 :                     , p[3]
     201           0 :                     , p[4]*cnv      };
     202             :   
     203           0 :   Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
     204           0 :   Double_t c32 =   x*cov[13] -     cov[ 8];
     205           0 :   Double_t c20 =   x*cov[10] -     cov[ 3];
     206           0 :   Double_t c21 =   x*cov[11] -     cov[ 4];
     207           0 :   Double_t c42 =   x*cov[14] -     cov[12];
     208             :   
     209           0 :   Double_t cc[15] = { cov[ 0]
     210           0 :                     , cov[ 1],     cov[ 2]
     211             :                     , c20,         c21,         c22
     212           0 :                     , cov[ 6],     cov[ 7],     c32,     cov[ 9]
     213           0 :                     , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
     214             :   
     215           0 :   Double_t mostProbablePt=AliExternalTrackParam::GetMostProbablePt();
     216           0 :   Double_t p0=TMath::Sign(1/mostProbablePt,pp[4]);
     217           0 :   Double_t w0=cc[14]/(cc[14] + p0*p0), w1=p0*p0/(cc[14] + p0*p0);
     218           0 :   AliDebug(3, Form("Pt mixing : w0[%4.2f] pt0[%5.3f] w1[%4.2f] pt[%5.3f]", w0, 1./p0, w1, 1./pp[4]));
     219             :   
     220           0 :   pp[4] = w0*p0 + w1*pp[4];
     221           0 :   cc[10]*=w1; cc[11]*=w1; cc[12]*=w1; cc[13]*=w1; cc[14]*=w1;
     222           0 :   Set(x,alpha,pp,cc);
     223           0 :   AliDebug(2, Form("Init @ x[%6.2f] pt[%5.3f]", x, 1./pp[4]));
     224             :   Int_t ncls = 0;
     225           0 :   for(int iplane=0; iplane<kNplane; iplane++){
     226           0 :     fTrackletIndex[iplane] = -1;
     227           0 :           if(!trklts[iplane].IsOK()) fTracklet[iplane] = NULL;
     228             :     else{ 
     229           0 :       fTracklet[iplane] = &trklts[iplane];
     230           0 :       ncls += fTracklet[iplane]->GetN();
     231             :     }
     232             :   }
     233           0 :   AliKalmanTrack::SetNumberOfClusters(ncls);            
     234           0 :   for(int i =0; i<3; i++) fBudget[i] = 0.;
     235             :   
     236             :   Float_t pid = 1./AliPID::kSPECIES;
     237           0 :   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
     238           0 :   SetLabel(-123456789); // reset label
     239           0 : }
     240             : 
     241             : //_______________________________________________________________
     242             : AliTRDtrackV1::~AliTRDtrackV1()
     243        1224 : {
     244             :   // Clean up all objects allocated by the track during its lifetime.
     245        1270 :   AliDebug(2, Form("Deleting track[%d]\n   fBackupTrack[%p] fTrackLow[%p] fTrackHigh[%p] Owner[%c].", fESDid, (void*)fBackupTrack, (void*)fTrackLow, (void*)fTrackHigh, TestBit(kOwner)?'y':'n'));
     246             : 
     247         508 :   if(fBackupTrack) delete fBackupTrack; fBackupTrack = NULL;
     248             : 
     249         676 :   if(fTrackLow) delete fTrackLow; fTrackLow = NULL;
     250         592 :   if(fTrackHigh) delete fTrackHigh; fTrackHigh = NULL;
     251             : 
     252        3556 :   for(Int_t ip=0; ip<kNplane; ip++){
     253        2560 :     if(TestBit(kOwner) && fTracklet[ip]) delete fTracklet[ip];
     254        1524 :     fTracklet[ip] = NULL;
     255        1524 :     fTrackletIndex[ip] = -1;
     256             :   }
     257         612 : }
     258             :         
     259             : //_______________________________________________________________
     260             : AliTRDtrackV1 &AliTRDtrackV1::operator=(const AliTRDtrackV1 &t)
     261             : {
     262             :   //
     263             :   // Assignment operator
     264             :   //
     265             : 
     266           0 :   if (this != &t) {
     267           0 :     AliKalmanTrack::operator=(t);
     268           0 :     ((AliTRDtrackV1 &) t).Copy(*this);
     269           0 :   }
     270             : 
     271           0 :   return *this;
     272             : 
     273             : }
     274             : 
     275             : //_____________________________________________________________________________
     276             : void AliTRDtrackV1::Copy(TObject &t) const
     277             : {
     278             :   //
     279             :   // Copy function
     280             :   //
     281             : 
     282           0 :   ((AliTRDtrackV1 &) t).fStatus         = fStatus;
     283           0 :   ((AliTRDtrackV1 &) t).fESDid          = fESDid;
     284           0 :   ((AliTRDtrackV1 &) t).fDE             = fDE;
     285           0 :   ((AliTRDtrackV1 &) t).fkReconstructor = fkReconstructor;
     286           0 :   ((AliTRDtrackV1 &) t).fBackupTrack    = 0x0;
     287           0 :   ((AliTRDtrackV1 &) t).fTrackLow       = 0x0;
     288           0 :   ((AliTRDtrackV1 &) t).fTrackHigh      = 0x0;
     289             : 
     290           0 :   for(Int_t ip = 0; ip < kNplane; ip++) { 
     291           0 :     ((AliTRDtrackV1 &) t).fTrackletIndex[ip] = fTrackletIndex[ip];
     292           0 :     ((AliTRDtrackV1 &) t).fTracklet[ip]      = fTracklet[ip];
     293             :   }
     294           0 :   if (fTrackLow) {
     295           0 :     ((AliTRDtrackV1 &) t).fTrackLow  = new AliExternalTrackParam(*fTrackLow);
     296           0 :   }
     297           0 :   if (fTrackHigh){
     298           0 :     ((AliTRDtrackV1 &) t).fTrackHigh = new AliExternalTrackParam(*fTrackHigh);
     299           0 :   }
     300             :  
     301           0 :   for (Int_t i = 0; i < 3; i++) {
     302           0 :     ((AliTRDtrackV1 &) t).fBudget[i] = fBudget[i];
     303             :   }
     304           0 :   for (Int_t is = 0; is < AliPID::kSPECIES; is++) {
     305           0 :     ((AliTRDtrackV1 &) t).fPID[is] = fPID[is];
     306             :   }  
     307             : 
     308           0 : }
     309             : 
     310             : //_______________________________________________________________
     311             : Int_t AliTRDtrackV1::CookLabel(Float_t wrong, Int_t *labs, Float_t *freq)
     312             : {
     313             : // Set MC label for this track
     314             : // On demand i.e. if arrays "labs" and "freq" are allocated by user returns :
     315             : //   nlabs = the no of distinct labels
     316             : //   labs  = array of distinct labels in decreasing order of frequency
     317             : //   freq  = frequency of each label  in decreasing order
     318             : 
     319             :   Int_t ncl(0);
     320         270 :   if(!(ncl = GetNumberOfClusters())) return 0;
     321             : 
     322          42 :   Int_t s[2][kMAXCLUSTERSPERTRACK];
     323       17724 :   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
     324        8820 :     s[0][i] = -1;
     325        8820 :     s[1][i] =  0;
     326             :   }
     327             : 
     328             :   Int_t label(-123456789), nlabels(0);
     329             :   AliTRDcluster *c(NULL);
     330         588 :   for (Int_t ip(0); ip < AliTRDgeometry::kNlayer; ip++) {
     331         458 :     if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
     332       25956 :     for (Int_t ic(0); ic < AliTRDseedV1::kNclusters; ic++) {
     333       12772 :       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
     334       35840 :       for (Int_t k(0); k < 3; k++) {
     335       13440 :         if ((label = c->GetLabel(k)) < 0) continue;
     336             :         Int_t j(0);
     337       92569 :         while(j < kMAXCLUSTERSPERTRACK){
     338      233345 :           if(s[0][j]!=label && s[1][j]!=0){j++; continue;}
     339        9680 :           if(!s[1][j]) nlabels++;
     340        9007 :           s[0][j] = label; s[1][j]++;
     341        9007 :           break;
     342             :         }
     343        9007 :       }
     344        4480 :     }
     345         206 :   }
     346             :   //printf("  Found %4d labels\n", nlabels);
     347             :   Float_t prob(1.);
     348          42 :   if(!nlabels){
     349           0 :     AliError(Form("No MC labels found for track %d.", fESDid));
     350           0 :     return 0;
     351          42 :   } else if(nlabels==1) {
     352           6 :     label = s[0][0];
     353           6 :     if(labs && freq){labs[0]=label; freq[0]=1.;}
     354             :   } else {
     355          36 :     Int_t idx[kMAXCLUSTERSPERTRACK];
     356          36 :     TMath::Sort(nlabels, s[1], idx);
     357          36 :     label = s[0][idx[0]]; prob = s[1][idx[0]]/Float_t(ncl);
     358          36 :     if(labs && freq){
     359           0 :       for (Int_t i(0); i<nlabels; i++){
     360           0 :         labs[i] = s[0][idx[i]];
     361           0 :         freq[i] = s[1][idx[i]]/Float_t(ncl);
     362             :       }
     363           0 :     }
     364          36 :   }
     365          42 :   SetLabel((1.-prob > wrong)?-label:label);
     366          42 :   return nlabels;
     367         146 : }
     368             : 
     369             : //_______________________________________________________________
     370             : Bool_t AliTRDtrackV1::CookPID()
     371             : {
     372             : //
     373             : // Cook the PID information for the track by delegating the omonim function of the tracklets. 
     374             : // Computes the number of tracklets used. The tracklet information are considered independent. 
     375             : // For the moment no global track measurement of PID is performed as for example to estimate 
     376             : // bremsstrahlung probability based on global chi2 of the track.
     377             : //
     378             : // The status bit AliESDtrack::kTRDpid is set during the call of AliTRDtrackV1::UpdateESDtrack().The PID performance of the 
     379             : //TRD for tracks with 6 tacklets is displayed below.
     380             : //Begin_Html
     381             : //<img src="TRD/trackPID.gif">
     382             : //End_Html
     383             : //
     384         208 :   const AliTRDPIDResponse *pidResponse = AliTRDcalibDB::Instance()->GetPIDResponse(fkReconstructor->GetRecoParam()->GetPIDmethod());
     385         104 :   if(!pidResponse){
     386           0 :     AliError("PID Response not available");
     387           0 :     return kFALSE;
     388             :   }
     389         104 :   Double_t dEdx[kNplane * (Int_t)AliTRDseedV1::kNdEdxSlices] = {0.};
     390         104 :   Float_t trackletP[kNplane] = {0.};
     391             : 
     392         104 :   fNdEdxSlices = pidResponse->GetNumberOfSlices();
     393        1456 :   for(Int_t iseed = 0; iseed < kNplane; iseed++){
     394         624 :     if(!fTracklet[iseed]) continue;
     395         206 :     trackletP[iseed] = fTracklet[iseed]->GetMomentum();
     396             : //    if(pidResponse->GetPIDmethod() == AliTRDPIDResponse::kLQ1D){
     397         206 :       dEdx[iseed] = fTracklet[iseed]->GetdQdl();
     398             : /*    } else {
     399             :       fTracklet[iseed]->CookdEdx(fNdEdxSlices);
     400             :       const Float_t *trackletdEdx = fTracklet[iseed]->GetdEdx();
     401             :       for(Int_t islice = 0; islice < fNdEdxSlices; islice++){
     402             :         dEdx[iseed*fNdEdxSlices + islice] = trackletdEdx[islice];
     403             :       }
     404             :     }*/
     405         206 :   }
     406         104 :   pidResponse->GetResponse(fNdEdxSlices, dEdx, trackletP, fPID);
     407             : 
     408             :   static Int_t nprint = 0;
     409         104 :   if(!nprint){
     410           2 :     AliTRDdEdxBaseUtils::PrintControl();
     411           2 :     nprint++;
     412           2 :   }
     413             : 
     414             :   //do truncated mean
     415         104 :   AliTRDdEdxCalibUtils::SetObjArray(AliTRDcalibDB::Instance()->GetPHQ());
     416         312 :   const Double_t mag    = AliTRDdEdxBaseUtils::IsExBOn() ? GetBz()  : -1;
     417         312 :   const Double_t charge = AliTRDdEdxBaseUtils::IsExBOn() ? Charge() : -1;
     418         104 :   fTruncatedMean = CookTruncatedMean(0, mag, charge, kTRUE, fNchamberdEdx, fNclusterdEdx);
     419             : 
     420             :   return kTRUE;
     421         208 : }
     422             : 
     423             : //___________________________________________________________
     424             : UChar_t AliTRDtrackV1::GetNumberOfTrackletsPID() const
     425             : {
     426             : // Retrieve number of tracklets used for PID calculation. 
     427             : 
     428             :   UChar_t nPID = 0;
     429        1560 :   for(int ip=0; ip<kNplane; ip++){
     430         830 :     if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
     431         206 :     if(!fTracklet[ip]->IsOK()) continue;
     432             :     
     433         206 :     nPID++;
     434         206 :   }
     435         104 :   return nPID;
     436             : }
     437             : 
     438             : //_______________________________________________________________
     439             : AliTRDcluster* AliTRDtrackV1::GetCluster(Int_t id)
     440             : {
     441             :   // Get the cluster at a certain position in the track
     442             :   Int_t n = 0;
     443           0 :   for(Int_t ip=0; ip<kNplane; ip++){
     444           0 :     if(!fTracklet[ip]) continue;
     445           0 :     if(n+fTracklet[ip]->GetN() <= id){ 
     446           0 :       n+=fTracklet[ip]->GetN();
     447           0 :       continue;
     448             :     }
     449             :     AliTRDcluster *c = NULL;
     450           0 :     for(Int_t ic=AliTRDseedV1::kNclusters; ic--;){
     451           0 :       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
     452             : 
     453           0 :       if(n<id){n++; continue;}
     454           0 :       return c;
     455             :     }
     456           0 :   }
     457           0 :   return NULL;
     458           0 : }
     459             : 
     460             : //_______________________________________________________________
     461             : Int_t  AliTRDtrackV1::GetClusterIndex(Int_t id) const
     462             : {
     463             :   // Get the cluster index at a certain position in the track
     464             :   Int_t n = 0;
     465           0 :   for(Int_t ip=0; ip<kNplane; ip++){
     466           0 :     if(!fTracklet[ip]) continue;
     467           0 :     if(n+fTracklet[ip]->GetN() <= id){ 
     468           0 :       n+=fTracklet[ip]->GetN();
     469           0 :       continue;
     470             :     }
     471           0 :     for(Int_t ic=AliTRDseedV1::kNclusters; ic--;){
     472           0 :       if(!(fTracklet[ip]->GetClusters(ic))) continue;
     473           0 :       if(n<id){n++; continue;}
     474           0 :       return fTracklet[ip]->GetIndexes(ic);
     475             :     }
     476             :   }
     477           0 :   return -1;
     478           0 : }
     479             : 
     480             : //_______________________________________________________________
     481             : Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt, Double_t *cov) const
     482             : {
     483             : // Compute chi2 between tracklet and track. The value is calculated at the radial position of the track
     484             : // equal to the reference radial position of the tracklet (see AliTRDseedV1)
     485             : // 
     486             : // The chi2 estimator is computed according to the following formula
     487             : // BEGIN_LATEX
     488             : // #chi^{2}=(X_{trklt}-X_{track})(C_{trklt}+C_{track})^{-1}(X_{trklt}-X_{track})^{T}
     489             : // END_LATEX
     490             : // where X=(y z), the position of the track/tracklet in the yz plane
     491             : // 
     492             : 
     493           0 :   Double_t p[2]   = { trklt->GetY(), trklt->GetZ()};
     494           0 :   trklt->GetCovAt(trklt->GetX(), cov);
     495           0 :   return AliExternalTrackParam::GetPredictedChi2(p, cov);
     496           0 : }
     497             : 
     498             : //_______________________________________________________________
     499             : Int_t AliTRDtrackV1::GetSector() const
     500             : {
     501        2402 :   return Int_t(GetAlpha()/AliTRDgeometry::GetAlpha() + (GetAlpha()>0. ? 0 : AliTRDgeometry::kNsector));
     502             : }
     503             : 
     504             : //_______________________________________________________________
     505             : Bool_t AliTRDtrackV1::IsEqual(const TObject *o) const
     506             : {
     507             :   // Checks whether two tracks are equal
     508         104 :   if (!o) return kFALSE;
     509         156 :   const AliTRDtrackV1 *inTrack = dynamic_cast<const AliTRDtrackV1*>(o);
     510         104 :   if (!inTrack) return kFALSE;
     511             :   
     512             :   //if ( fPIDquality != inTrack->GetPIDquality() ) return kFALSE;
     513             :   
     514           0 :   if(memcmp(fPID, inTrack->fPID, AliPID::kSPECIES*sizeof(Double32_t))) return kFALSE;
     515           0 :   if(memcmp(fBudget, inTrack->fBudget, 3*sizeof(Double32_t))) return kFALSE;
     516           0 :   if(memcmp(&fDE, &inTrack->fDE, sizeof(Double32_t))) return kFALSE;
     517           0 :   if(memcmp(&fFakeRatio, &inTrack->fFakeRatio, sizeof(Double32_t))) return kFALSE;
     518           0 :   if(memcmp(&fChi2, &inTrack->fChi2, sizeof(Double32_t))) return kFALSE;
     519           0 :   if(memcmp(&fMass, &inTrack->fMass, sizeof(Double32_t))) return kFALSE;
     520           0 :   if( fLab != inTrack->fLab ) return kFALSE;
     521           0 :   if( fN != inTrack->fN ) return kFALSE;
     522           0 :   Double32_t l(0.), in(0.);
     523           0 :   l = GetIntegratedLength(); in = inTrack->GetIntegratedLength();
     524           0 :   if(memcmp(&l, &in, sizeof(Double32_t))) return kFALSE;
     525           0 :   l=GetX(); in=inTrack->GetX();
     526           0 :   if(memcmp(&l, &in, sizeof(Double32_t))) return kFALSE;
     527           0 :   l = GetAlpha(); in = inTrack->GetAlpha();
     528           0 :   if(memcmp(&l, &in, sizeof(Double32_t))) return kFALSE;
     529           0 :   if(memcmp(GetParameter(), inTrack->GetParameter(), 5*sizeof(Double32_t))) return kFALSE;
     530           0 :   if(memcmp(GetCovariance(), inTrack->GetCovariance(), 15*sizeof(Double32_t))) return kFALSE;
     531             :   
     532           0 :   for (Int_t iTracklet = 0; iTracklet < kNplane; iTracklet++){
     533           0 :     AliTRDseedV1 *curTracklet = fTracklet[iTracklet];
     534           0 :     AliTRDseedV1 *inTracklet = inTrack->GetTracklet(iTracklet);
     535           0 :     if (curTracklet && inTracklet){
     536           0 :       if (! curTracklet->IsEqual(inTracklet) ) {
     537           0 :         curTracklet->Print();
     538           0 :         inTracklet->Print();
     539           0 :         return kFALSE;
     540             :       }
     541             :     } else {
     542             :       // if one tracklet exists, and corresponding 
     543             :       // in other track doesn't - return kFALSE
     544           0 :       if(inTracklet || curTracklet) return kFALSE;
     545             :     }
     546           0 :   }
     547             : 
     548           0 :   return kTRUE;
     549          52 : }
     550             : 
     551             : //_______________________________________________________________
     552             : Bool_t AliTRDtrackV1::IsElectron() const
     553             : {
     554           0 :   if(GetPID(0) > fkReconstructor->GetRecoParam()->GetPIDThreshold(GetP())) return kTRUE;
     555           0 :   return kFALSE;
     556           0 : }
     557             : 
     558             :         
     559             : //_____________________________________________________________________________
     560             : Int_t AliTRDtrackV1::MakeBackupTrack()
     561             : {
     562             : //
     563             : // Creates a backup track
     564             : // TO DO update quality check of the track.
     565             : //
     566             : 
     567             :   Float_t occupancy(0.); Int_t n(0), ncls(0);
     568        2538 :   for(Int_t il(AliTRDgeometry::kNlayer); il--;){ 
     569        1236 :     if(!fTracklet[il]) continue;
     570         684 :     n++; 
     571         684 :     occupancy+=fTracklet[il]->GetTBoccupancy()/AliTRDseedV1::kNtb;
     572         684 :     ncls += fTracklet[il]->GetN();
     573             :   }
     574         206 :   if(!n) return -1;
     575         206 :   occupancy/=n;
     576             : 
     577             :   //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);  
     578             :   
     579             :   Int_t failedCutId(0);
     580         217 :   if(GetChi2()/n > 5.0) failedCutId=1; 
     581         412 :   if(occupancy < 0.7) failedCutId=2;
     582             :   //if(ratio1 >   0.6) && 
     583             :   //if(ratio0+ratio1  >   1.5) && 
     584         262 :   if(GetNCross() != 0)  failedCutId=3;
     585         206 :   if(TMath::Abs(GetSnp()) > 0.85) failedCutId=4;
     586         217 :   if(ncls < 20) failedCutId=5;
     587             : 
     588         206 :   if(failedCutId){ 
     589         618 :     AliDebug(2, Form("\n"
     590             :       "chi2/tracklet < 5.0   [%c] %5.2f\n"
     591             :       "occupancy > 0.7       [%c] %4.2f\n"
     592             :       "NCross == 0           [%c] %d\n"
     593             :       "Abs(snp) < 0.85       [%c] %4.2f\n"
     594             :       "NClusters > 20        [%c] %d"
     595             :       ,(GetChi2()/n<5.0)?'y':'n', GetChi2()/n
     596             :       ,(occupancy>0.7)?'y':'n', occupancy
     597             :       ,(GetNCross()==0)?'y':'n', GetNCross()
     598             :       ,(TMath::Abs(GetSnp())<0.85)?'y':'n', TMath::Abs(GetSnp())
     599             :       ,(ncls>20)?'y':'n', ncls
     600             :     ));
     601         206 :     return failedCutId;
     602             :   }
     603             : 
     604           0 :   if(fBackupTrack) {
     605           0 :     fBackupTrack->~AliTRDtrackV1();
     606           0 :     new(fBackupTrack) AliTRDtrackV1((AliTRDtrackV1&)(*this));
     607           0 :     return 0;
     608             :   }
     609           0 :   fBackupTrack = new AliTRDtrackV1((AliTRDtrackV1&)(*this));
     610           0 :   return 0;
     611         206 : }
     612             : 
     613             : //_____________________________________________________________________________
     614             : Int_t AliTRDtrackV1::GetProlongation(Double_t xk, Double_t &y, Double_t &z) const
     615             : {
     616             :   //
     617             :   // Find a prolongation at given x
     618             :   // Return -1 if it does not exist
     619             :   //  
     620             : 
     621       23030 :   Double_t bz = GetBz();
     622       11541 :   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return -1;
     623       11489 :   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return -1;
     624             : 
     625       11489 :   return 1;  
     626             : 
     627       11515 : }
     628             : 
     629             : //_____________________________________________________________________________
     630             : Bool_t AliTRDtrackV1::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
     631             : {
     632             :   //
     633             :   // Propagates this track to a reference plane defined by "xk" [cm] 
     634             :   // correcting for the mean crossed material.
     635             :   //
     636             :   // "xx0"  - thickness/rad.length [units of the radiation length] 
     637             :   // "xrho" - thickness*density    [g/cm^2] 
     638             :   // 
     639             : 
     640       20278 :   if (TMath::Abs(xk - GetX())<AliTRDReconstructor::GetEpsilon()*0.1) return kTRUE; // 10% of the tracker precision
     641             : 
     642       10139 :   Double_t xyz0[3] = {GetX(), GetY(), GetZ()}, // track position BEFORE propagation 
     643             :            b[3];    // magnetic field 
     644       10139 :   GetBxByBz(b);
     645       10139 :   if(!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE;
     646             :  
     647             :   // local track position AFTER propagation 
     648       10139 :   Double_t xyz1[3] = {GetX(), GetY(), GetZ()};
     649             : //  printf("x0[%6.2f] -> x1[%6.2f] dx[%6.2f] rho[%f]\n", xyz0[0], xyz1[0], xyz0[0]-xk, xrho/TMath::Abs(xyz0[0]-xk));
     650       10139 :   if(xyz0[0] < xk) {
     651        7307 :     xrho = -xrho;
     652        7307 :     if (IsStartedTimeIntegral()) {
     653       14614 :       Double_t l2  = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) 
     654        7307 :                                + (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) 
     655        7307 :                                + (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
     656        7307 :       Double_t crv = AliExternalTrackParam::GetC(b[2]);
     657        7307 :       if (TMath::Abs(l2*crv) > 0.0001) {
     658             :         // Make correction for curvature if neccesary
     659        6933 :         l2 = 0.5 * TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) 
     660             :                              + (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]));
     661        6933 :         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
     662        6933 :         l2 = TMath::Sqrt(l2*l2 + (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
     663        6933 :       }
     664        7307 :       AddTimeStep(l2);
     665        7307 :     }
     666             :   }
     667       10142 :   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0, xrho, fMass)) return kFALSE;
     668             : 
     669             : 
     670             :   {
     671             : 
     672             :     // Energy losses
     673       10136 :     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
     674       10136 :     if (fMass<0) p2 *= 4; // q=2
     675       10136 :     Double_t beta2 = p2 / (p2 + fMass*fMass);
     676       20272 :     if ((beta2 < 1.0e-10) || 
     677       10136 :         ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
     678           0 :       return kFALSE;
     679             :     }
     680             : 
     681       10136 :     Double_t dE    = 0.153e-3 / beta2 
     682       10136 :                    * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
     683       10136 :                    * xrho;
     684       10136 :     if (fMass<0) dE *= 4; // q=2
     685       10136 :     fBudget[0] += xrho;
     686             : 
     687             :     /*
     688             :     // Suspicious part - think about it ?
     689             :     Double_t kinE =  TMath::Sqrt(p2);
     690             :     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //      
     691             :     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch 
     692             :     */
     693             :  
     694       10136 :     fDE += dE;
     695             : 
     696             :     /*
     697             :     // Suspicious ! I.B.
     698             :     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation 
     699             :     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+fMass*fMass)/(p2*p2);
     700             :     fCcc += sigmac2;
     701             :     fCee += fX*fX * sigmac2;  
     702             :     */
     703             : 
     704       10136 :   }
     705             : 
     706       10136 :   return kTRUE;
     707       20278 : }
     708             : 
     709             : //_____________________________________________________________________________
     710             : Int_t   AliTRDtrackV1::PropagateToR(Double_t r,Double_t step)
     711             : {
     712             :   //
     713             :   // Propagate track to the radial position
     714             :   // Rotation always connected to the last track position
     715             :   //
     716             : 
     717           0 :   Double_t xyz0[3];
     718           0 :   Double_t xyz1[3];
     719           0 :   Double_t y;
     720           0 :   Double_t z; 
     721             : 
     722           0 :   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
     723             :   // Direction +-
     724           0 :   Double_t dir    = (radius > r) ? -1.0 : 1.0;   
     725             : 
     726           0 :   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
     727             : 
     728           0 :     GetXYZ(xyz0);       
     729           0 :     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
     730           0 :     Rotate(alpha,kTRUE);
     731           0 :     GetXYZ(xyz0);       
     732           0 :     if(GetProlongation(x,y,z)<0) return -1;
     733           0 :     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
     734           0 :     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
     735           0 :     xyz1[2] = z;
     736           0 :     Double_t param[7];
     737           0 :     if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param)<=0.) return -1;
     738           0 :     if (param[1] <= 0) {
     739           0 :       param[1] = 100000000;
     740           0 :     }
     741           0 :     PropagateTo(x,param[1],param[0]*param[4]);
     742             : 
     743           0 :   } 
     744             : 
     745           0 :   GetXYZ(xyz0); 
     746           0 :   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
     747           0 :   Rotate(alpha,kTRUE);
     748           0 :   GetXYZ(xyz0); 
     749           0 :   if(GetProlongation(r,y,z)<0) return -1;
     750           0 :   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
     751           0 :   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
     752           0 :   xyz1[2] = z;
     753           0 :   Double_t param[7];
     754           0 :   if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param) <= 0.) return -1;
     755             : 
     756           0 :   if (param[1] <= 0) {
     757           0 :     param[1] = 100000000;
     758           0 :   }
     759           0 :   PropagateTo(r,param[1],param[0]*param[4]);
     760             : 
     761           0 :   return 0;
     762             : 
     763           0 : }
     764             : 
     765             : //_____________________________________________________________________________
     766             : void AliTRDtrackV1::Print(Option_t *o) const
     767             : {
     768             :   // Print track status
     769           0 :   AliInfo(Form("PID [%4.1f %4.1f %4.1f %4.1f %4.1f]", 1.E2*fPID[0], 1.E2*fPID[1], 1.E2*fPID[2], 1.E2*fPID[3], 1.E2*fPID[4]));
     770           0 :   AliInfo(Form("Material[%5.2f %5.2f %5.2f]", fBudget[0], fBudget[1], fBudget[2]));
     771             : 
     772           0 :   AliInfo(Form("x[%7.2f] t[%7.4f] alpha[%f] mass[%f]", GetX(), GetIntegratedLength(), GetAlpha(), fMass));
     773           0 :   AliInfo(Form("Ntr[%1d] NtrPID[%1d] Ncl[%3d] lab[%3d]", GetNumberOfTracklets(), GetNumberOfTrackletsPID(), fN, fLab));
     774             : 
     775           0 :   printf("|X| = (");
     776           0 :   const Double_t *curP = GetParameter();
     777           0 :   for (Int_t i = 0; i < 5; i++) printf("%7.2f ", curP[i]);
     778           0 :   printf(")\n");
     779             : 
     780           0 :   printf("|V| = \n");
     781           0 :   const Double_t *curC = GetCovariance();
     782           0 :   for (Int_t i = 0, j=4, k=0; i<15; i++, k++){
     783           0 :     printf("%7.2f ", curC[i]);
     784           0 :     if(k==j){ 
     785           0 :       printf("\n");
     786           0 :       k=-1; j--;
     787           0 :     }
     788             :   }
     789           0 :   if(strcmp(o, "a")!=0) return;
     790             : 
     791           0 :   for(Int_t ip=0; ip<kNplane; ip++){
     792           0 :     if(!fTracklet[ip]) continue;
     793           0 :     fTracklet[ip]->Print(o);
     794           0 :   }
     795           0 : }
     796             : 
     797             : 
     798             : //_____________________________________________________________________________
     799             : Bool_t AliTRDtrackV1::Rotate(Double_t alpha, Bool_t absolute)
     800             : {
     801             :   //
     802             :   // Rotates track parameters in R*phi plane
     803             :   // if absolute rotation alpha is in global system
     804             :   // otherwise alpha rotation is relative to the current rotation angle
     805             :   //  
     806             : 
     807          46 :   if (absolute) alpha -= GetAlpha();
     808             :   //else fNRotate++;
     809             : 
     810          46 :   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
     811             : }
     812             : 
     813             : //___________________________________________________________
     814             : void AliTRDtrackV1::SetNumberOfClusters() 
     815             : {
     816             : // Calculate the number of clusters attached to this track
     817             :         
     818             :   Int_t ncls = 0;
     819        6180 :   for(int ip=0; ip<kNplane; ip++){
     820        5208 :     if(fTracklet[ip] && fTrackletIndex[ip] >= 0) ncls += fTracklet[ip]->GetN();
     821             :   }
     822         412 :   AliKalmanTrack::SetNumberOfClusters(ncls);    
     823         412 : }
     824             : 
     825             :         
     826             : //_______________________________________________________________
     827             : void AliTRDtrackV1::SetOwner()
     828             : {
     829             :   //
     830             :   // Toggle ownership of tracklets
     831             :   //
     832             : 
     833         208 :   if(TestBit(kOwner)) return;
     834        1456 :   for (Int_t ip = 0; ip < kNplane; ip++) {
     835         830 :     if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
     836         412 :     fTracklet[ip] = new AliTRDseedV1(*fTracklet[ip]);
     837         206 :     fTracklet[ip]->SetOwner();
     838         206 :   }
     839         104 :   SetBit(kOwner);
     840         208 : }
     841             : 
     842             : //_______________________________________________________________
     843             : void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *const trklt, Int_t index)
     844             : {
     845             :   //
     846             :   // Set the tracklets
     847             :   //
     848         824 :   Int_t plane = trklt->GetPlane();
     849             : 
     850         412 :   fTracklet[plane]      = trklt;
     851         412 :   fTrackletIndex[plane] = index;
     852         412 : }
     853             : 
     854             : //_______________________________________________________________
     855             : void AliTRDtrackV1::SetTrackIn()
     856             : {
     857             : //  Save location of birth for the TRD track
     858             : //  If the pointer is not valid allocate memory
     859             : //
     860          84 :   const AliExternalTrackParam *op = dynamic_cast<const AliExternalTrackParam*>(this);
     861             : 
     862             :   //printf("SetTrackIn() : fTrackLow[%p]\n", (void*)fTrackLow);
     863          42 :   if(fTrackLow){
     864           0 :     fTrackLow->~AliExternalTrackParam();
     865           0 :     new(fTrackLow) AliExternalTrackParam(*op);
     866          84 :   } else fTrackLow = new AliExternalTrackParam(*op);
     867          42 : }
     868             : 
     869             : //_______________________________________________________________
     870             : void AliTRDtrackV1::SetTrackOut(const AliExternalTrackParam *op)
     871             : {
     872             : //  Save location of death for the TRD track
     873             : //  If the pointer is not valid allocate memory
     874             : //
     875          84 :   if(!op) op = dynamic_cast<const AliExternalTrackParam*>(this);
     876          42 :   if(fTrackHigh){
     877           0 :     fTrackHigh->~AliExternalTrackParam();
     878           0 :     new(fTrackHigh) AliExternalTrackParam(*op);
     879          84 :   } else fTrackHigh = new AliExternalTrackParam(*op);
     880          42 : }
     881             : 
     882             : //_______________________________________________________________
     883             : void AliTRDtrackV1::UnsetTracklet(Int_t plane)
     884             : {
     885           0 :   if(plane<0) return;
     886           0 :   fTrackletIndex[plane] = -1;
     887           0 :   fTracklet[plane] = NULL;
     888           0 : }
     889             : 
     890             : 
     891             : //_______________________________________________________________
     892             : void AliTRDtrackV1::UpdateChi2(Float_t chi2)
     893             : {
     894             : // Update chi2/track with one tracklet contribution
     895         824 :   SetChi2(GetChi2() + chi2);
     896         412 : }
     897             : 
     898             : //_______________________________________________________________
     899             : void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
     900             : {
     901             :   //
     902             :   // Update the TRD PID information in the ESD track
     903             :   //
     904             : 
     905             : //   Int_t nslices = AliTRDcalibDB::Instance()->GetPIDResponse(fkReconstructor->GetRecoParam()->GetPIDmethod())->GetNumberOfSlices();
     906             :   // number of tracklets used for PID calculation
     907         208 :   UChar_t nPID = GetNumberOfTrackletsPID();
     908             :   // number of tracklets attached to the track
     909         104 :   UChar_t nTrk = GetNumberOfTracklets();
     910             :   // pack the two numbers together and store them in the ESD
     911         104 :   track->SetTRDntracklets(nPID | (nTrk<<3));
     912             :   // allocate space to store raw PID signals dEdx & momentum
     913             :   // independent of the method used to calculate PID (see below AliTRDPIDResponse)
     914         104 :   track->SetNumberOfTRDslices((AliTRDseedV1::kNdEdxSlices+2)*AliTRDgeometry::kNlayer);
     915             :   // store raw signals
     916         104 :   Float_t p, sp; Double_t spd;
     917        1456 :   for (Int_t ip = 0; ip < kNplane; ip++) {
     918         830 :     if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
     919             :     // Fill TRD dEdx info into ESD track
     920             :     // a. Set Summed dEdx into the first slice
     921         206 :     track->SetTRDslice(fTracklet[ip]->GetdQdl(), ip, 0); 
     922             :     // b. Set NN dEdx slices
     923         206 :     fTracklet[ip]->CookdEdx(AliTRDPIDResponse::kNslicesNN);
     924         206 :     const Float_t *dedx = fTracklet[ip]->GetdEdx();
     925        3502 :     for (Int_t js(0), ks(1); js < AliTRDPIDResponse::kNslicesNN; js++, ks++, dedx++){
     926        1442 :       if(ks>=AliTRDseedV1::kNdEdxSlices){
     927           0 :         AliError(Form("Exceed allocated space for dEdx slices."));
     928           0 :         break;
     929             :       }
     930        1442 :       track->SetTRDslice(*dedx, ip, ks);
     931             :     }
     932             :     // fill TRD momentum info into ESD track
     933         206 :     p = fTracklet[ip]->GetMomentum(&sp); spd = sp;
     934         206 :     track->SetTRDmomentum(p, ip, &spd);
     935             :     // store global quality per tracklet instead of momentum error
     936             :     // 26.09.11 A.Bercuci
     937             :     // first implementation store no. of time bins filled in tracklet (5bits  see "y" bits) and
     938             :     // no. of double clusters in case of pad row cross (4bits see "x" bits)
     939             :     // bit map for tracklet quality xxxxyyyyy
     940             :     // 27.10.11 A.Bercuci
     941             :     // add chamber status bit "z" bit
     942             :     // bit map for tracklet quality zxxxxyyyyy
     943             :     // 12.11.11 A.Bercuci
     944             :     // fit tracklet quality into the field fTRDTimeBin [Char_t]
     945             :     // bit map for tracklet quality zxxyyyyy
     946             :     // The information should be retrieved by the following functions of AliESDtrack for each TRD layer
     947             :     // GetTRDtrkltOccupancy(layer) -> no of TB filled in tracklet
     948             :     // GetTRDtrkltClCross(layer)   -> no of TB filled in crossing pad rows
     949             :     // IsTRDtrkltChmbGood(layer)   -> status of the chamber from which the tracklet is found
     950         426 :     Int_t nCross(fTracklet[ip]->IsRowCross()?fTracklet[ip]->GetTBcross():0); if(nCross>3) nCross = 3;
     951         206 :     Char_t trackletQ = Char_t(fTracklet[ip]->GetTBoccupancy() | (nCross<<5) | (fTracklet[ip]->IsChmbGood()<<7));
     952         206 :     track->SetTRDTimBin(trackletQ, ip);
     953         206 :   }
     954             :   // store PID probabilities
     955         104 :   track->SetTRDpid(fPID);
     956             : 
     957             :   //store truncated mean
     958         104 :   track->SetTRDsignal(fTruncatedMean);
     959         104 :   track->SetTRDNchamberdEdx(fNchamberdEdx);
     960         104 :   track->SetTRDNclusterdEdx(fNclusterdEdx);
     961         104 : }
     962             : 
     963             : //_______________________________________________________________
     964             : Double_t  AliTRDtrackV1::CookTruncatedMean(const Bool_t kinvq, const Double_t mag, const Int_t charge, const Int_t kcalib, Int_t &nch, Int_t &ncls, TVectorD *Qs, TVectorD *Xs, Int_t timeBin0, Int_t timeBin1, Int_t tstep) const
     965             : {
     966             :   //
     967             :   //Origin: Xianguo Lu <xianguo.lu@cern.ch>, Marian Ivanov <marian.ivanov@cern.ch>
     968             :   //
     969             : 
     970         104 :   TVectorD arrayQ(200), arrayX(200);
     971         208 :   ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, this, timeBin0, timeBin1, tstep);
     972             : 
     973         312 :   const TObjArray *cobj = kcalib ? AliTRDdEdxCalibUtils::GetObj(kinvq, mag, charge) : NULL;
     974             :   
     975         104 :   const Double_t tmean = AliTRDdEdxReconUtils::ToyCook(kinvq, ncls, &arrayQ, &arrayX, cobj);
     976             : 
     977         208 :   nch = AliTRDdEdxReconUtils::UpdateArrayX(ncls, &arrayX);
     978             : 
     979         104 :   if(Qs && Xs){
     980           0 :     (*Qs)=arrayQ;
     981           0 :     (*Xs)=arrayX;
     982             :   }
     983             : 
     984             :   //printf("\ntest %.10f %d %d\n", tmean, nch, ncls);
     985             : 
     986             :   return tmean;
     987         104 : }
     988             : 
     989             : //_______________________________________________________________
     990             : TObject* AliTRDtrackV1::Clone(const char* newname) const
     991             : {
     992             :   // temporary override TObject::Clone to avoid crashes in reco
     993             :   AliTRDtrackV1* src = (AliTRDtrackV1*)this;
     994           0 :   Bool_t isown = src->IsOwner();
     995             :   //AliInfo(Form("src_owner %d",isown));
     996           0 :   AliTRDtrackV1* dst = new AliTRDtrackV1(*src);
     997           0 :   if (isown) {
     998           0 :     src->SetBit(kOwner);
     999           0 :     dst->SetOwner();
    1000           0 :   }
    1001           0 :   return dst;
    1002           0 : }

Generated by: LCOV version 1.11