TRD/AliTRDseedV1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  The TRD track seed                                                    //
  21 //                                                                        //
  22 //  Authors:                                                              //
  23 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
  24 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  25 //                                                                        //
  26 ////////////////////////////////////////////////////////////////////////////
  27
  28 #include "TMath.h"
  29 #include "TLinearFitter.h"
  30
  31 #include "AliLog.h"
  32 #include "AliMathBase.h"
  33
  34 #include "AliTRDseedV1.h"
  35 #include "AliTRDcluster.h"
  36 #include "AliTRDtrack.h"
  37 #include "AliTRDcalibDB.h"
  38 #include "AliTRDstackLayer.h"
  39 #include "AliTRDrecoParam.h"
  40 #include "AliTRDgeometry.h"
  41 #include "Cal/AliTRDCalPID.h"
  42
  43 #define SEED_DEBUG
  44
  45 ClassImp(AliTRDseedV1)
  46
  47 //____________________________________________________________________
  48 AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p)
  49   :AliTRDseed()
  50   ,fPlane(layer)
  51   ,fOwner(kFALSE)
  52   ,fMom(0.)
  53   ,fRecoParam(p)
  54 {
  55   //
  56   // Constructor
  57   //
  58         for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
  59         for(int itb=0; itb < knTimebins; itb++){
  60                 fdQdl[itb]  = 0.;
  61                 fdQ[itb]    = 0.;
  62         }
  63         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
  64 }
  65
  66 //____________________________________________________________________
  67 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
  68   :AliTRDseed((AliTRDseed&)ref)
  69   ,fPlane(ref.fPlane)
  70   ,fOwner(kFALSE)
  71   ,fMom(ref.fMom)
  72   ,fRecoParam(ref.fRecoParam)
  73 {
  74   //
  75   // Copy Constructor performing a deep copy
  76   //
  77
  78         //AliInfo("");
  79         if(ref.fOwner) SetOwner();
  80         for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
  81         for(int itb=0; itb < knTimebins; itb++){
  82                 fdQdl[itb] = ref.fdQdl[itb];
  83                 fdQ[itb]   = ref.fdQ[itb];
  84         }
  85         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
  86 }
  87
  88
  89 //____________________________________________________________________
  90 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
  91 {
  92   //
  93   // Assignment Operator using the copy function
  94   //
  95
  96         //AliInfo("");
  97         if(this != &ref){
  98                 ref.Copy(*this);
  99         }
 100         return *this;
 101
 102 }
 103
 104 //____________________________________________________________________
 105 AliTRDseedV1::~AliTRDseedV1()
 106 {
 107   //
 108   // Destructor. The RecoParam object belongs to the underlying tracker.
 109   //
 110
 111         //AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO"));
 112
 113         if(fOwner)
 114                 for(int itb=0; itb<knTimebins; itb++){
 115                         if(!fClusters[itb]) continue;
 116                         //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
 117                         delete fClusters[itb];
 118                         fClusters[itb] = 0x0;
 119                 }
 120 }
 121
 122 //____________________________________________________________________
 123 void AliTRDseedV1::Copy(TObject &ref) const
 124 {
 125   //
 126   // Copy function
 127   //
 128
 129         //AliInfo("");
 130         AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
 131
 132         target.fPlane         = fPlane;
 133         target.fMom           = fMom;
 134         target.fRecoParam     = fRecoParam;
 135
 136         for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
 137         for(int itb=0; itb < knTimebins; itb++){
 138                 target.fdQdl[itb] = fdQdl[itb];
 139                 target.fdQ[itb]   = fdQ[itb];
 140         }
 141         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
 142
 143         AliTRDseed::Copy(target);
 144 }
 145
 146
 147 //____________________________________________________________
 148 void AliTRDseedV1::Init(AliTRDtrack *track)
 149 {
 150 // Initialize this tracklet using the track information
 151 //
 152 // Parameters:
 153 //   track - the TRD track used to initialize the tracklet
 154 //
 155 // Detailed description
 156 // The function sets the starting point and direction of the
 157 // tracklet according to the information from the TRD track.
 158 //
 159 // Caution
 160 // The TRD track has to be propagated to the beginning of the
 161 // chamber where the tracklet will be constructed
 162 //
 163
 164         Double_t y, z;
 165         track->GetProlongation(fX0, y, z);
 166         fYref[0] = y;
 167         fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
 168         fZref[0] = z;
 169         fZref[1] = track->GetTgl();
 170
 171         //printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl());
 172 }
 173
 174 //____________________________________________________________________
 175 Double_t* AliTRDseedV1::GetProbability()
 176 {
 177 // Fill probability array for tracklet from the DB.
 178 //
 179 // Parameters
 180 //
 181 // Output
 182 //   returns pointer to the probability array and 0x0 if missing DB access
 183 //
 184 // Detailed description
 185
 186
 187         // retrive calibration db
 188   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 189   if (!calibration) {
 190     AliError("No access to calibration data");
 191     return 0x0;
 192   }
 193
 194   // Retrieve the CDB container class with the parametric detector response
 195   const AliTRDCalPID *pd = calibration->GetPIDObject(fRecoParam->GetPIDMethod());
 196   if (!pd) {
 197     AliError("No access to AliTRDCalPID object");
 198     return 0x0;
 199   }
 200
 201         // calculate tracklet length TO DO
 202   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 203   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
 204
 205   //calculate dE/dx
 206   CookdEdx(fRecoParam->GetNdEdxSlices());
 207
 208   // Sets the a priori probabilities
 209   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 210     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, fPlane);
 211   }
 212
 213         return &fProb[0];
 214 }
 215
 216 //____________________________________________________________________
 217 void AliTRDseedV1::CookdEdx(Int_t nslices)
 218 {
 219 // Calculates average dE/dx for all slices and store them in the internal array fdEdx.
 220 //
 221 // Parameters:
 222 //  nslices : number of slices for which dE/dx should be calculated
 223 // Output:
 224 //
 225 // Detailed description
 226 // Calculates average dE/dx for all slices. Depending on the PID methode
 227 // the number of slices can be 3 (LQ) or 8(NN). The calculation is based
 228 // on previously calculated quantities dQ/dl of each cluster. The
 229 // following effects are included in the calculation:
 230 // 1. calibration values for t0 and vdrift
 231 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
 232 //
 233
 234         Int_t nclusters[knSlices];
 235         for(int i=0; i<knSlices; i++){
 236                 fdEdx[i]     = 0.;
 237                 nclusters[i] = 0;
 238         }
 239
 240         AliTRDcluster *cluster = 0x0;
 241         for(int ic=0; ic<fgTimeBins; ic++){
 242                 if(!(cluster = fClusters[ic])) continue;
 243                 Int_t tb = cluster->GetLocalTimeBin();
 244
 245                 // consider calibration effects
 246                 if(tb < fTimeBin0 || tb >= fTimeBin0+fTimeBinsRange) continue;
 247
 248                 // consider cluster sharing ... TO DO
 249                 //if(fRecoParam->GetClusterSharing() && cluster->GetSharing()) continue;
 250
 251                 Int_t slice = (tb-fTimeBin0)*nslices/fTimeBinsRange;
 252                 fdEdx[slice]   += fdQdl[ic];
 253                 nclusters[slice]++;
 254         } // End of loop over clusters
 255
 256         // calculate mean charge per slice
 257         for(int is=0; is<nslices; is++) if(nclusters[is]) fdEdx[is] /= nclusters[is];
 258 }
 259
 260 //____________________________________________________________________
 261 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
 262 {
 263   //
 264   // Returns a quality measurement of the current seed
 265   //
 266
 267         Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 268         return .5 * (18.0 - fN2)
 269                 + 10.* TMath::Abs(fYfit[1] - fYref[1])
 270                 + 5.* TMath::Abs(fYfit[0] - fYref[0] + zcorr)
 271                 + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
 272 }
 273
 274 //____________________________________________________________________
 275 void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
 276 {
 277 // Computes covariance in the y-z plane at radial point x
 278
 279         const Float_t k0= .2; // to be checked in FindClusters
 280         Double_t sy20   = k0*TMath::Tan(fYfit[1]); sy20 *= sy20;
 281
 282         Double_t sy2    = fSigmaY2*fSigmaY2 + sy20;
 283         Double_t sz2    = fPadLength/12.;
 284
 285         //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
 286
 287         cov[0] = sy2;
 288         cov[1] = fTilt*(sy2-sz2);
 289         cov[2] = sz2;
 290 }
 291
 292 //____________________________________________________________________
 293 void AliTRDseedV1::SetdQdl(Double_t length)
 294 {
 295         for(int ic=0; ic<fgTimeBins; ic++) fdQdl[ic] = fdQ[ic] *length;
 296 }
 297
 298 //____________________________________________________________________
 299 void AliTRDseedV1::SetOwner(Bool_t own)
 300 {
 301         //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
 302
 303         if(own){
 304                 for(int ic=0; ic<knTimebins; ic++){
 305                         if(!fClusters[ic]) continue;
 306                         fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
 307                 }
 308                 fOwner = kTRUE;
 309         } else {
 310                 if(fOwner){
 311                         for(int ic=0; ic<knTimebins; ic++){
 312                                 if(!fClusters[ic]) continue;
 313                                 delete fClusters[ic];
 314                                 //fClusters[ic] = tracker->GetClusters(index) TODO
 315                         }
 316                 }
 317                 fOwner = kFALSE;
 318         }
 319 }
 320
 321 //____________________________________________________________________
 322 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
 323                                        , Float_t quality
 324                                        , Bool_t kZcorr
 325                                        , AliTRDcluster *c)
 326 {
 327   //
 328   // Iterative process to register clusters to the seed.
 329   // In iteration 0 we try only one pad-row and if quality not
 330   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
 331   //
 332
 333         if(!fRecoParam){
 334                 AliError("Seed can not be used without a valid RecoParam.");
 335                 return kFALSE;
 336         }
 337
 338         //AliInfo(Form("TimeBins = %d TimeBinsRange = %d", fgTimeBins, fTimeBinsRange));
 339
 340         Float_t  tquality;
 341         Double_t kroady = fRecoParam->GetRoad1y();
 342         Double_t kroadz = fPadLength * .5 + 1.;
 343
 344         // initialize configuration parameters
 345         Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 346         Int_t   niter = kZcorr ? 1 : 2;
 347
 348         Double_t yexp, zexp;
 349         Int_t ncl = 0;
 350         // start seed update
 351         for (Int_t iter = 0; iter < niter; iter++) {
 352         //AliInfo(Form("iter = %i", iter));
 353                 ncl = 0;
 354                 for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
 355                         // define searching configuration
 356                         Double_t dxlayer = layer[iTime].GetX() - fX0;
 357                         if(c){
 358                                 zexp = c->GetZ();
 359                                 //Try 2 pad-rows in second iteration
 360                                 if (iter > 0) {
 361                                         zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
 362                                         if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
 363                                         if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
 364                                 }
 365                         } else zexp = fZref[0];
 366                         yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
 367                         // get  cluster
 368 //                      printf("xexp = %3.3f ,yexp = %3.3f, zexp = %3.3f\n",layer[iTime].GetX(),yexp,zexp);
 369 //                      printf("layer[%i].GetNClusters() = %i\n", iTime, layer[iTime].GetNClusters());
 370                         Int_t    index = layer[iTime].SearchNearestCluster(yexp, zexp, kroady, kroadz);
 371
 372 //                      printf("%d[%d] x[%7.3f | %7.3f] y[%7.3f] z[%7.3f]\n", iTime, layer[iTime].GetNClusters(), dxlayer, layer[iTime].GetX(), yexp, zexp);
 373 //                      for(Int_t iclk = 0; iclk < layer[iTime].GetNClusters(); iclk++){
 374 //                              AliTRDcluster *testcl = layer[iTime].GetCluster(iclk);
 375 //                              printf("Cluster %i: %d x = %7.3f, y = %7.3f, z = %7.3f\n", iclk, testcl->GetLocalTimeBin(), testcl->GetX(), testcl->GetY(), testcl->GetZ());
 376 //                      }
 377 //                      printf("Index = %i\n",index);
 378
 379                         if (index < 0) continue;
 380
 381                         // Register cluster
 382                         AliTRDcluster *cl = (AliTRDcluster*) layer[iTime].GetCluster(index);
 383
 384                         //printf("Cluster %i(0x%x): x = %3.3f, y = %3.3f, z = %3.3f\n", index, cl, cl->GetX(), cl->GetY(), cl->GetZ());
 385
 386                         Int_t globalIndex = layer[iTime].GetGlobalIndex(index);
 387                         fIndexes[iTime]  = globalIndex;
 388                         fClusters[iTime] = cl;
 389                         fX[iTime]        = dxlayer;
 390                         fY[iTime]        = cl->GetY();
 391                         fZ[iTime]        = cl->GetZ();
 392                         fdQ[iTime]       = cl->GetQ()/layer[iTime].GetdX();
 393
 394                         // Debugging
 395                         ncl++;
 396                 }
 397
 398 #ifdef SEED_DEBUG
 399 //              Int_t nclusters = 0;
 400 //              Float_t fD[iter] = 0.;
 401 //              for(int ic=0; ic<fgTimeBins+1; ic++){
 402 //                      AliTRDcluster *ci = fClusters[ic];
 403 //                      if(!ci) continue;
 404 //                      for(int jc=ic+1; jc<fgTimeBins+1; jc++){
 405 //                              AliTRDcluster *cj = fClusters[jc];
 406 //                              if(!cj) continue;
 407 //                              fD[iter] += TMath::Sqrt((ci->GetY()-cj->GetY())*(ci->GetY()-cj->GetY())+
 408 //                              (ci->GetZ()-cj->GetZ())*(ci->GetZ()-cj->GetZ()));
 409 //                              nclusters++;
 410 //                      }
 411 //              }
 412 //              if(nclusters) fD[iter] /= float(nclusters);
 413 #endif
 414
 415                 AliTRDseed::Update();
 416
 417                 if(IsOK()){
 418                         tquality = GetQuality(kZcorr);
 419                         if(tquality < quality) break;
 420                         else quality = tquality;
 421                 }
 422                 kroadz *= 2.;
 423         } // Loop: iter
 424         if (!IsOK()) return kFALSE;
 425
 426         CookLabels();
 427         UpdateUsed();
 428         return kTRUE;
 429 }
 430
 431 //____________________________________________________________________
 432 Bool_t  AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
 433                                        ,Bool_t kZcorr)
 434 {
 435   //
 436   // Projective algorithm to attach clusters to seeding tracklets
 437   //
 438   // Parameters
 439   //
 440   // Output
 441   //
 442   // Detailed description
 443   // 1. Collapse x coordinate for the full detector plane
 444   // 2. truncated mean on y (r-phi) direction
 445   // 3. purge clusters
 446   // 4. truncated mean on z direction
 447   // 5. purge clusters
 448   // 6. fit tracklet
 449   //
 450
 451         if(!fRecoParam){
 452                 AliError("Seed can not be used without a valid RecoParam.");
 453                 return kFALSE;
 454         }
 455
 456         const Int_t kClusterCandidates = 2 * knTimebins;
 457
 458         //define roads
 459         Double_t kroady = fRecoParam->GetRoad1y();
 460         Double_t kroadz = fPadLength * 1.5 + 1.;
 461         // correction to y for the tilting angle
 462         Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 463
 464         // working variables
 465         AliTRDcluster *clusters[kClusterCandidates];
 466         Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
 467                 yres[kClusterCandidates], zres[kClusterCandidates];
 468         Int_t ncl, *index = 0x0, tboundary[knTimebins];
 469
 470         // Do cluster projection
 471         Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
 472         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
 473                 fX[iTime] = layer[iTime].GetX() - fX0;
 474                 zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
 475                 yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
 476
 477                 // build condition and process clusters
 478                 cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
 479                 cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
 480                 layer[iTime].GetClusters(cond, index, ncl);
 481                 for(Int_t ic = 0; ic<ncl; ic++){
 482                         AliTRDcluster *c = layer[iTime].GetCluster(index[ic]);
 483                         clusters[nYclusters] = c;
 484                         yres[nYclusters++] = c->GetY() - yexp[iTime];
 485                         if(nYclusters >= kClusterCandidates) {
 486                                 AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
 487                                 kEXIT = kTRUE;
 488                                 break;
 489                         }
 490                 }
 491                 tboundary[iTime] = nYclusters;
 492                 if(kEXIT) break;
 493         }
 494
 495         // Evaluate truncated mean on the y direction
 496         Double_t mean, sigma;
 497         AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
 498         //purge cluster candidates
 499         Int_t nZclusters = 0;
 500         for(Int_t ic = 0; ic<nYclusters; ic++){
 501                 if(yres[ic] - mean > 4. * sigma){
 502                         clusters[ic] = 0x0;
 503                         continue;
 504                 }
 505                 zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
 506         }
 507
 508         // Evaluate truncated mean on the z direction
 509         AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
 510         //purge cluster candidates
 511         for(Int_t ic = 0; ic<nZclusters; ic++){
 512                 if(zres[ic] - mean > 4. * sigma){
 513                         clusters[ic] = 0x0;
 514                         continue;
 515                 }
 516         }
 517
 518
 519         // Select only one cluster/TimeBin
 520         Int_t lastCluster = 0;
 521         fN2 = 0;
 522         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
 523                 ncl = tboundary[iTime] - lastCluster;
 524                 if(!ncl) continue;
 525                 AliTRDcluster *c = 0x0;
 526                 if(ncl == 1){
 527                         c = clusters[lastCluster];
 528                 } else if(ncl > 1){
 529                         Float_t dold = 9999.; Int_t iptr = lastCluster;
 530                         for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
 531                                 if(!clusters[ic]) continue;
 532                                 Float_t y = yexp[iTime] - clusters[ic]->GetY();
 533                                 Float_t z = zexp[iTime] - clusters[ic]->GetZ();
 534                                 Float_t d = y * y + z * z;
 535                                 if(d > dold) continue;
 536                                 dold = d;
 537                                 iptr = ic;
 538                         }
 539                         c = clusters[iptr];
 540                 }
 541                 //Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index);
 542                 //fIndexes[iTime]  = GlobalIndex;
 543                 fClusters[iTime] = c;
 544                 fY[iTime]        = c->GetY();
 545                 fZ[iTime]        = c->GetZ();
 546                 fdQ[iTime]       = c->GetQ()/layer[iTime].GetdX();
 547                 lastCluster      = tboundary[iTime];
 548                 fN2++;
 549         }
 550
 551         // number of minimum numbers of clusters expected for the tracklet
 552         Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fgTimeBins);
 553   if (fN2 < kClmin){
 554                 AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
 555     fN2 = 0;
 556     return kFALSE;
 557   }
 558
 559         // update used clusters
 560         fNUsed = 0;
 561         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
 562                 if(!fClusters[iTime]) continue;
 563                 if((fClusters[iTime]->IsUsed())) fNUsed++;
 564         }
 565
 566   if (fN2-fNUsed < kClmin){
 567                 AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
 568     fN2 = 0;
 569     return kFALSE;
 570   }
 571
 572         return kTRUE;
 573 }
 574
 575 //____________________________________________________________________
 576 Bool_t AliTRDseedV1::Fit()
 577 {
 578   //
 579   // Linear fit of the tracklet
 580   //
 581   // Parameters :
 582   //
 583   // Output :
 584   //  True if successful
 585   //
 586   // Detailed description
 587   // 2. Check if tracklet crosses pad row boundary
 588   // 1. Calculate residuals in the y (r-phi) direction
 589   // 3. Do a Least Square Fit to the data
 590   //
 591
 592         //Float_t  sigmaexp  = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
 593   Float_t  ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
 594   Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
 595
 596         // calculate residuals
 597         Float_t yres[knTimebins]; // y (r-phi) residuals
 598         Int_t zint[knTimebins],   // Histograming of the z coordinate
 599               zout[2*knTimebins];//
 600
 601         fN = 0;
 602         for (Int_t iTime = 0; iTime < fTimeBinsRange; iTime++) {
 603     if (!fClusters[iTime]) continue;
 604     yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime];
 605                 zint[fN] = Int_t(fZ[iTime]);
 606                 fN++;
 607         }
 608
 609         // calculate pad row boundary crosses
 610         Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBinsRange);
 611         Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE);
 612   fZProb   = zout[0];
 613   if(nz <= 1) zout[3] = 0;
 614   if(zout[1] + zout[3] < kClmin) {
 615                 AliWarning(Form("Not enough clusters to fit the cross boundary tracklet %d [%d].", zout[1]+zout[3], kClmin));
 616                 return kFALSE;
 617         }
 618   // Z distance bigger than pad - length
 619   if (TMath::Abs(zout[0]-zout[2]) > fPadLength) zout[3]=0;
 620
 621
 622   Double_t sumw   = 0.,
 623                 sumwx  = 0.,
 624                 sumwx2 = 0.,
 625                 sumwy  = 0.,
 626                 sumwxy = 0.,
 627                 sumwz  = 0.,
 628                 sumwxz = 0.;
 629         Int_t npads;
 630   fMPads = 0;
 631         fMeanz = 0.;
 632         // we will use only the clusters which are in the detector range
 633         for(int iTime=0; iTime<fTimeBinsRange; iTime++){
 634     fUsable[iTime] = kFALSE;
 635     if (!fClusters[iTime]) continue;
 636                 npads = fClusters[iTime]->GetNPads();
 637
 638                 fUsable[iTime] = kTRUE;
 639     fN2++;
 640     fMPads += npads;
 641     Float_t weight = 1.0;
 642     if(npads > 5) weight = 0.2;
 643     else if(npads > 4) weight = 0.5;
 644     sumw   += weight;
 645     sumwx  += fX[iTime] * weight;
 646     sumwx2 += fX[iTime] * fX[iTime] * weight;
 647     sumwy  += weight * yres[iTime];
 648     sumwxy += weight * yres[iTime] * fX[iTime];
 649     sumwz  += weight * fZ[iTime];
 650     sumwxz += weight * fZ[iTime] * fX[iTime];
 651         }
 652   if (fN2 < kClmin){
 653                 AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
 654     fN2 = 0;
 655     return kFALSE;
 656   }
 657   fMeanz = sumwz / sumw;
 658         fNChange = 0;
 659
 660         // Tracklet on boundary
 661   Float_t correction = 0;
 662   if (fNChange > 0) {
 663     if (fMeanz < fZProb) correction =  ycrosscor;
 664     if (fMeanz > fZProb) correction = -ycrosscor;
 665   }
 666
 667   Double_t det = sumw * sumwx2 - sumwx * sumwx;
 668   fYfitR[0]    = (sumwx2 * sumwy  - sumwx * sumwxy) / det;
 669   fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
 670
 671   fSigmaY2 = 0;
 672   for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 673     if (!fUsable[i]) continue;
 674     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
 675     fSigmaY2 += delta*delta;
 676   }
 677   fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
 678
 679   fZfitR[0]  = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
 680   fZfitR[1]  = (sumw   * sumwxz - sumwx * sumwz)  / det;
 681   fZfit[0]   = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
 682   fZfit[1]   = (sumw   * sumwxz - sumwx * sumwz)  / det;
 683   fYfitR[0] += fYref[0] + correction;
 684   fYfitR[1] += fYref[1];
 685   fYfit[0]   = fYfitR[0];
 686   fYfit[1]   = fYfitR[1];
 687
 688         return kTRUE;
 689 }
 690
 691 //_____________________________________________________________________________
 692 Float_t AliTRDseedV1::FitRiemanTilt(AliTRDseedV1 *cseed, Bool_t terror)
 693 {
 694   //
 695   // Fit the Rieman tilt
 696   //
 697
 698   // Fitting with tilting pads - kz not fixed
 699         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 700         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 701   TLinearFitter fitterT2(4,"hyp4");
 702   fitterT2.StoreData(kTRUE);
 703   Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
 704
 705   Int_t npointsT = 0;
 706   fitterT2.ClearPoints();
 707
 708   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 709 //              printf("\nLayer %d\n", iLayer);
 710 //     cseed[iLayer].Print();
 711                 if (!cseed[iLayer].IsOK()) continue;
 712     Double_t tilt = cseed[iLayer].fTilt;
 713
 714     for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
 715 //                      printf("\ttime %d\n", itime);
 716       if (!cseed[iLayer].fUsable[itime]) continue;
 717       // x relative to the midle chamber
 718       Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
 719       Double_t y = cseed[iLayer].fY[itime];
 720       Double_t z = cseed[iLayer].fZ[itime];
 721
 722       //
 723       // Tilted rieman
 724       //
 725       Double_t uvt[6];
 726       Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0;      // Global x
 727       Double_t t  = 1.0 / (x2*x2 + y*y);
 728       uvt[1]  = t;
 729       uvt[0]  = 2.0 * x2   * uvt[1];
 730       uvt[2]  = 2.0 * tilt * uvt[1];
 731       uvt[3]  = 2.0 * tilt *uvt[1] * x;
 732       uvt[4]  = 2.0 * (y + tilt * z) * uvt[1];
 733
 734       Double_t error = 2.0 * uvt[1];
 735       if (terror) {
 736         error *= cseed[iLayer].fSigmaY;
 737       }
 738       else {
 739         error *= 0.2; //Default error
 740       }
 741 //                      printf("\tadd point :\n");
 742 //                      for(int i=0; i<5; i++) printf("%f ", uvt[i]);
 743 //                      printf("\n");
 744       fitterT2.AddPoint(uvt,uvt[4],error);
 745       npointsT++;
 746
 747     }
 748
 749   }
 750   fitterT2.Eval();
 751   Double_t rpolz0 = fitterT2.GetParameter(3);
 752   Double_t rpolz1 = fitterT2.GetParameter(4);
 753
 754   //
 755   // Linear fitter  - not possible to make boundaries
 756   // non accept non possible z and dzdx combination
 757   //
 758   Bool_t acceptablez = kTRUE;
 759   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 760     if (cseed[iLayer].IsOK()) {
 761       Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
 762       if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
 763         acceptablez = kFALSE;
 764       }
 765     }
 766   }
 767   if (!acceptablez) {
 768     Double_t zmf  = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
 769     Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
 770     fitterT2.FixParameter(3,zmf);
 771     fitterT2.FixParameter(4,dzmf);
 772     fitterT2.Eval();
 773     fitterT2.ReleaseParameter(3);
 774     fitterT2.ReleaseParameter(4);
 775     rpolz0 = fitterT2.GetParameter(3);
 776     rpolz1 = fitterT2.GetParameter(4);
 777   }
 778
 779   Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
 780   Double_t params[3];
 781   params[0] =  fitterT2.GetParameter(0);
 782   params[1] =  fitterT2.GetParameter(1);
 783   params[2] =  fitterT2.GetParameter(2);
 784   Double_t curvature =  1.0 + params[1] * params[1] - params[2] * params[0];
 785
 786   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 787
 788     Double_t  x  = cseed[iLayer].fX0;
 789     Double_t  y  = 0;
 790     Double_t  dy = 0;
 791     Double_t  z  = 0;
 792     Double_t  dz = 0;
 793
 794     // y
 795     Double_t res2 = (x * params[0] + params[1]);
 796     res2 *= res2;
 797     res2  = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
 798     if (res2 >= 0) {
 799       res2 = TMath::Sqrt(res2);
 800       y    = (1.0 - res2) / params[0];
 801     }
 802
 803     //dy
 804     Double_t x0 = -params[1] / params[0];
 805     if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
 806       Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
 807       if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
 808         Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
 809         if (params[0] < 0) res *= -1.0;
 810         dy = res;
 811       }
 812     }
 813     z  = rpolz0 + rpolz1 * (x - xref2);
 814     dz = rpolz1;
 815     cseed[iLayer].fYref[0] = y;
 816     cseed[iLayer].fYref[1] = dy;
 817     cseed[iLayer].fZref[0] = z;
 818     cseed[iLayer].fZref[1] = dz;
 819     cseed[iLayer].fC       = curvature;
 820
 821   }
 822
 823   return chi2TR;
 824
 825 }
 826
 827 //___________________________________________________________________
 828 void AliTRDseedV1::Print()
 829 {
 830   //
 831   // Printing the seedstatus
 832   //
 833
 834         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 835         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 836
 837         printf("Seed status :\n");
 838         printf("  fTilt      = %f\n", fTilt);
 839         printf("  fPadLength = %f\n", fPadLength);
 840         printf("  fX0        = %f\n", fX0);
 841         for(int ic=0; ic<nTimeBins; ic++) {
 842           const Char_t *isUsable = fUsable[ic]?"Yes":"No";
 843           printf("  %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n"
 844                 , ic
 845                 , fX[ic]
 846                 , fY[ic]
 847                 , fZ[ic]
 848                 , fIndexes[ic]
 849                 , ((void*) fClusters[ic])
 850                 , isUsable);
 851         }
 852
 853         printf("  fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
 854         printf("  fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
 855         printf("  fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
 856         printf("  fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]);
 857         printf("  fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
 858         printf("  fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]);
 859         printf("  fSigmaY =%f\n", fSigmaY);
 860         printf("  fSigmaY2=%f\n", fSigmaY2);
 861         printf("  fMeanz  =%f\n", fMeanz);
 862         printf("  fZProb  =%f\n", fZProb);
 863         printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
 864         printf("  fN      =%d\n", fN);
 865         printf("  fN2     =%d (>8 isOK)\n",fN2);
 866         printf("  fNUsed  =%d\n", fNUsed);
 867         printf("  fFreq   =%d\n", fFreq);
 868         printf("  fNChange=%d\n",  fNChange);
 869         printf("  fMPads  =%f\n", fMPads);
 870
 871         printf("  fC      =%f\n", fC);
 872         printf("  fCC     =%f\n",fCC);
 873         printf("  fChi2   =%f\n", fChi2);
 874         printf("  fChi2Z  =%f\n", fChi2Z);
 875
 876 }