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