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 #include "TClonesArray.h" // tmp
  31 #include <TTreeStream.h>
  32
  33 #include "AliLog.h"
  34 #include "AliMathBase.h"
  35
  36 #include "AliTRDpadPlane.h"
  37 #include "AliTRDcluster.h"
  38 #include "AliTRDseedV1.h"
  39 #include "AliTRDtrackV1.h"
  40 #include "AliTRDcalibDB.h"
  41 #include "AliTRDchamberTimeBin.h"
  42 #include "AliTRDtrackingChamber.h"
  43 #include "AliTRDtrackerV1.h"
  44 #include "AliTRDReconstructor.h"
  45 #include "AliTRDrecoParam.h"
  46 #include "Cal/AliTRDCalPID.h"
  47
  48 ClassImp(AliTRDseedV1)
  49
  50 //____________________________________________________________________
  51 AliTRDseedV1::AliTRDseedV1(Int_t det)
  52   :AliTRDseed()
  53   ,fReconstructor(0x0)
  54   ,fClusterIter(0x0)
  55   ,fClusterIdx(0)
  56   ,fDet(det)
  57   ,fMom(0.)
  58   ,fSnp(0.)
  59   ,fTgl(0.)
  60   ,fdX(0.)
  61   ,fXref(0.)
  62 {
  63   //
  64   // Constructor
  65   //
  66   //printf("AliTRDseedV1::AliTRDseedV1()\n");
  67
  68   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
  69   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
  70   fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
  71 }
  72
  73 //____________________________________________________________________
  74 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
  75   :AliTRDseed((AliTRDseed&)ref)
  76   ,fReconstructor(ref.fReconstructor)
  77   ,fClusterIter(0x0)
  78   ,fClusterIdx(0)
  79   ,fDet(ref.fDet)
  80   ,fMom(ref.fMom)
  81   ,fSnp(ref.fSnp)
  82   ,fTgl(ref.fTgl)
  83   ,fdX(ref.fdX)
  84   ,fXref(ref.fXref)
  85 {
  86   //
  87   // Copy Constructor performing a deep copy
  88   //
  89
  90   //printf("AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &)\n");
  91   SetBit(kOwner, kFALSE);
  92   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
  93   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
  94   memcpy(fRefCov, ref.fRefCov, 3*sizeof(Double_t));
  95 }
  96
  97
  98 //____________________________________________________________________
  99 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
 100 {
 101   //
 102   // Assignment Operator using the copy function
 103   //
 104
 105   if(this != &ref){
 106     ref.Copy(*this);
 107   }
 108   SetBit(kOwner, kFALSE);
 109
 110   return *this;
 111
 112 }
 113
 114 //____________________________________________________________________
 115 AliTRDseedV1::~AliTRDseedV1()
 116 {
 117   //
 118   // Destructor. The RecoParam object belongs to the underlying tracker.
 119   //
 120
 121   //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
 122
 123   if(IsOwner())
 124     for(int itb=0; itb<knTimebins; itb++){
 125       if(!fClusters[itb]) continue;
 126       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
 127       delete fClusters[itb];
 128       fClusters[itb] = 0x0;
 129     }
 130 }
 131
 132 //____________________________________________________________________
 133 void AliTRDseedV1::Copy(TObject &ref) const
 134 {
 135   //
 136   // Copy function
 137   //
 138
 139   //AliInfo("");
 140   AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
 141
 142   target.fClusterIter   = 0x0;
 143   target.fClusterIdx    = 0;
 144   target.fDet           = fDet;
 145   target.fMom           = fMom;
 146   target.fSnp           = fSnp;
 147   target.fTgl           = fTgl;
 148   target.fdX            = fdX;
 149   target.fXref          = fXref;
 150   target.fReconstructor = fReconstructor;
 151
 152   for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
 153   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
 154   memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
 155
 156   AliTRDseed::Copy(target);
 157 }
 158
 159
 160 //____________________________________________________________
 161 Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
 162 {
 163 // Initialize this tracklet using the track information
 164 //
 165 // Parameters:
 166 //   track - the TRD track used to initialize the tracklet
 167 //
 168 // Detailed description
 169 // The function sets the starting point and direction of the
 170 // tracklet according to the information from the TRD track.
 171 //
 172 // Caution
 173 // The TRD track has to be propagated to the beginning of the
 174 // chamber where the tracklet will be constructed
 175 //
 176
 177   Double_t y, z;
 178   if(!track->GetProlongation(fX0, y, z)) return kFALSE;
 179   fYref[0] = y;
 180   fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
 181   fZref[0] = z;
 182   fZref[1] = track->GetTgl();
 183
 184   const Double_t *cov = track->GetCovariance();
 185   fRefCov[0] = cov[0]; // Var(y)
 186   fRefCov[1] = cov[1]; // Cov(yz)
 187   fRefCov[2] = cov[5]; // Var(z)
 188
 189   //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());
 190   return kTRUE;
 191 }
 192
 193
 194 //____________________________________________________________________
 195 void AliTRDseedV1::CookdEdx(Int_t nslices)
 196 {
 197 // Calculates average dE/dx for all slices and store them in the internal array fdEdx.
 198 //
 199 // Parameters:
 200 //  nslices : number of slices for which dE/dx should be calculated
 201 // Output:
 202 //  store results in the internal array fdEdx. This can be accessed with the method
 203 //  AliTRDseedV1::GetdEdx()
 204 //
 205 // Detailed description
 206 // Calculates average dE/dx for all slices. Depending on the PID methode
 207 // the number of slices can be 3 (LQ) or 8(NN).
 208 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t)) i.e.
 209 //
 210 // dQ/dl = qc/(dx * sqrt(1 + dy/dx^2 + dz/dx^2))
 211 //
 212 // The following effects are included in the calculation:
 213 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
 214 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
 215 // 3. cluster size
 216 //
 217
 218   Int_t nclusters[knSlices];
 219   for(int i=0; i<knSlices; i++){
 220     fdEdx[i]     = 0.;
 221     nclusters[i] = 0;
 222   }
 223   Float_t clength = (/*.5 * */AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 224
 225   AliTRDcluster *cluster = 0x0;
 226   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
 227     if(!(cluster = fClusters[ic])) continue;
 228     Float_t x = cluster->GetX();
 229
 230     // Filter clusters for dE/dx calculation
 231
 232     // 1.consider calibration effects for slice determination
 233     Int_t slice;
 234     if(cluster->IsInChamber()) slice = Int_t(TMath::Abs(fX0 - x) * nslices / clength);
 235     else slice = x < fX0 ? 0 : nslices-1;
 236
 237     // 2. take sharing into account
 238     Float_t w = cluster->IsShared() ? .5 : 1.;
 239
 240     // 3. take into account large clusters TODO
 241     //w *= c->GetNPads() > 3 ? .8 : 1.;
 242
 243     //CHECK !!!
 244     fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
 245     nclusters[slice]++;
 246   } // End of loop over clusters
 247
 248   //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
 249   if(nslices == AliTRDReconstructor::kLQslices){
 250   // calculate mean charge per slice (only LQ PID)
 251     for(int is=0; is<nslices; is++){
 252       if(nclusters[is]) fdEdx[is] /= nclusters[is];
 253     }
 254   }
 255 }
 256
 257
 258 //____________________________________________________________________
 259 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
 260 {
 261   return fClusters[ic] ? TMath::Abs(fClusters[ic]->GetQ()) /fdX / TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]) : 0.;
 262 }
 263
 264 //____________________________________________________________________
 265 Double_t* AliTRDseedV1::GetProbability()
 266 {
 267 // Fill probability array for tracklet from the DB.
 268 //
 269 // Parameters
 270 //
 271 // Output
 272 //   returns pointer to the probability array and 0x0 if missing DB access
 273 //
 274 // Detailed description
 275
 276
 277   // retrive calibration db
 278   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 279   if (!calibration) {
 280     AliError("No access to calibration data");
 281     return 0x0;
 282   }
 283
 284   if (!fReconstructor) {
 285     AliError("Reconstructor not set.");
 286     return 0x0;
 287   }
 288
 289   // Retrieve the CDB container class with the parametric detector response
 290   const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
 291   if (!pd) {
 292     AliError("No access to AliTRDCalPID object");
 293     return 0x0;
 294   }
 295   //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
 296
 297   // calculate tracklet length TO DO
 298   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 299   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
 300
 301   //calculate dE/dx
 302   CookdEdx(fReconstructor->GetNdEdxSlices());
 303
 304   // Sets the a priori probabilities
 305   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 306     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());
 307   }
 308
 309   return &fProb[0];
 310 }
 311
 312 //____________________________________________________________________
 313 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
 314 {
 315   //
 316   // Returns a quality measurement of the current seed
 317   //
 318
 319   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 320   return
 321       .5 * TMath::Abs(18.0 - fN2)
 322     + 10.* TMath::Abs(fYfit[1] - fYref[1])
 323     + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
 324     + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
 325 }
 326
 327 //____________________________________________________________________
 328 void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
 329 {
 330 // Computes covariance in the y-z plane at radial point x
 331
 332   Int_t ic = 0; while (!fClusters[ic]) ic++;
 333   AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
 334   Double_t exB         = fCalib->GetOmegaTau(fCalib->GetVdriftAverage(fClusters[ic]->GetDetector()), -AliTracker::GetBz()*0.1);
 335
 336   Double_t sy2    = fSigmaY2*fSigmaY2 + .2*(fYfit[1]-exB)*(fYfit[1]-exB);
 337   Double_t sz2    = fPadLength/12.;
 338
 339
 340   //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
 341
 342   cov[0] = sy2;
 343   cov[1] = fTilt*(sy2-sz2);
 344   cov[2] = sz2;
 345
 346   // insert systematic uncertainties calibration and misalignment
 347   Double_t sys[15];
 348   fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
 349   cov[0] += (sys[0]*sys[0]);
 350   cov[2] += (sys[1]*sys[1]);
 351 }
 352
 353
 354 //____________________________________________________________________
 355 void AliTRDseedV1::SetOwner()
 356 {
 357   //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
 358
 359   if(TestBit(kOwner)) return;
 360   for(int ic=0; ic<knTimebins; ic++){
 361     if(!fClusters[ic]) continue;
 362     fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
 363   }
 364   SetBit(kOwner);
 365 }
 366
 367 //____________________________________________________________________
 368 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
 369 {
 370   //
 371   // Iterative process to register clusters to the seed.
 372   // In iteration 0 we try only one pad-row and if quality not
 373   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
 374   //
 375   // debug level 7
 376   //
 377
 378   if(!fReconstructor->GetRecoParam() ){
 379     AliError("Seed can not be used without a valid RecoParam.");
 380     return kFALSE;
 381   }
 382
 383   AliTRDchamberTimeBin *layer = 0x0;
 384   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
 385     AliTRDtrackingChamber ch(*chamber);
 386     ch.SetOwner();
 387     (*AliTRDtrackerV1::DebugStreamer()) << "AttachClustersIter"
 388       << "chamber.="   << &ch
 389       << "tracklet.="  << this
 390       << "\n";
 391   }
 392
 393   Float_t  tquality;
 394   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
 395   Double_t kroadz = fPadLength * .5 + 1.;
 396
 397   // initialize configuration parameters
 398   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 399   Int_t   niter = kZcorr ? 1 : 2;
 400
 401   Double_t yexp, zexp;
 402   Int_t ncl = 0;
 403   // start seed update
 404   for (Int_t iter = 0; iter < niter; iter++) {
 405     ncl = 0;
 406     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 407       if(!(layer = chamber->GetTB(iTime))) continue;
 408       if(!Int_t(*layer)) continue;
 409
 410       // define searching configuration
 411       Double_t dxlayer = layer->GetX() - fX0;
 412       if(c){
 413         zexp = c->GetZ();
 414         //Try 2 pad-rows in second iteration
 415         if (iter > 0) {
 416           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
 417           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
 418           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
 419         }
 420       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
 421       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
 422
 423       // Get and register cluster
 424       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
 425       if (index < 0) continue;
 426       AliTRDcluster *cl = (*layer)[index];
 427
 428       fIndexes[iTime]  = layer->GetGlobalIndex(index);
 429       fClusters[iTime] = cl;
 430       fY[iTime]        = cl->GetY();
 431       fZ[iTime]        = cl->GetZ();
 432       ncl++;
 433     }
 434     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
 435
 436     if(ncl>1){
 437       // calculate length of the time bin (calibration aware)
 438       Int_t irp = 0; Float_t x[2]; Int_t tb[2];
 439       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 440         if(!fClusters[iTime]) continue;
 441         x[irp]  = fClusters[iTime]->GetX();
 442         tb[irp] = iTime;
 443         irp++;
 444         if(irp==2) break;
 445       }
 446       fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
 447
 448       // update X0 from the clusters (calibration/alignment aware)
 449       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 450         if(!(layer = chamber->GetTB(iTime))) continue;
 451         if(!layer->IsT0()) continue;
 452         if(fClusters[iTime]){
 453           fX0 = fClusters[iTime]->GetX();
 454           break;
 455         } else { // we have to infere the position of the anode wire from the other clusters
 456           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
 457             if(!fClusters[jTime]) continue;
 458             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
 459             break;
 460           }
 461         }
 462       }
 463
 464       // update YZ reference point
 465       // TODO
 466
 467       // update x reference positions (calibration/alignment aware)
 468       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 469         if(!fClusters[iTime]) continue;
 470         fX[iTime] = fClusters[iTime]->GetX() - fX0;
 471       }
 472
 473       AliTRDseed::Update();
 474     }
 475     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
 476
 477     if(IsOK()){
 478       tquality = GetQuality(kZcorr);
 479       if(tquality < quality) break;
 480       else quality = tquality;
 481     }
 482     kroadz *= 2.;
 483   } // Loop: iter
 484   if (!IsOK()) return kFALSE;
 485
 486   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
 487   UpdateUsed();
 488   return kTRUE;
 489 }
 490
 491 //____________________________________________________________________
 492 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
 493                                       ,Bool_t kZcorr)
 494 {
 495   //
 496   // Projective algorithm to attach clusters to seeding tracklets
 497   //
 498   // Parameters
 499   //
 500   // Output
 501   //
 502   // Detailed description
 503   // 1. Collapse x coordinate for the full detector plane
 504   // 2. truncated mean on y (r-phi) direction
 505   // 3. purge clusters
 506   // 4. truncated mean on z direction
 507   // 5. purge clusters
 508   // 6. fit tracklet
 509   //
 510
 511   if(!fReconstructor->GetRecoParam() ){
 512     AliError("Seed can not be used without a valid RecoParam.");
 513     return kFALSE;
 514   }
 515
 516   const Int_t kClusterCandidates = 2 * knTimebins;
 517
 518   //define roads
 519   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
 520   Double_t kroadz = fPadLength * 1.5 + 1.;
 521   // correction to y for the tilting angle
 522   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 523
 524   // working variables
 525   AliTRDcluster *clusters[kClusterCandidates];
 526   Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
 527     yres[kClusterCandidates], zres[kClusterCandidates];
 528   Int_t ncl, *index = 0x0, tboundary[knTimebins];
 529
 530   // Do cluster projection
 531   AliTRDchamberTimeBin *layer = 0x0;
 532   Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
 533   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 534     if(!(layer = chamber->GetTB(iTime))) continue;
 535     if(!Int_t(*layer)) continue;
 536
 537     fX[iTime] = layer->GetX() - fX0;
 538     zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
 539     yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
 540
 541     // build condition and process clusters
 542     cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
 543     cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
 544     layer->GetClusters(cond, index, ncl);
 545     for(Int_t ic = 0; ic<ncl; ic++){
 546       AliTRDcluster *c = layer->GetCluster(index[ic]);
 547       clusters[nYclusters] = c;
 548       yres[nYclusters++] = c->GetY() - yexp[iTime];
 549       if(nYclusters >= kClusterCandidates) {
 550         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
 551         kEXIT = kTRUE;
 552         break;
 553       }
 554     }
 555     tboundary[iTime] = nYclusters;
 556     if(kEXIT) break;
 557   }
 558
 559   // Evaluate truncated mean on the y direction
 560   Double_t mean, sigma;
 561   AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
 562   // purge cluster candidates
 563   Int_t nZclusters = 0;
 564   for(Int_t ic = 0; ic<nYclusters; ic++){
 565     if(yres[ic] - mean > 4. * sigma){
 566       clusters[ic] = 0x0;
 567       continue;
 568     }
 569     zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
 570   }
 571
 572   // Evaluate truncated mean on the z direction
 573   AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
 574   // purge cluster candidates
 575   for(Int_t ic = 0; ic<nZclusters; ic++){
 576     if(zres[ic] - mean > 4. * sigma){
 577       clusters[ic] = 0x0;
 578       continue;
 579     }
 580   }
 581
 582
 583   // Select only one cluster/TimeBin
 584   Int_t lastCluster = 0;
 585   fN2 = 0;
 586   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 587     ncl = tboundary[iTime] - lastCluster;
 588     if(!ncl) continue;
 589     Int_t iptr = lastCluster;
 590     if(ncl > 1){
 591       Float_t dold = 9999.;
 592       for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
 593         if(!clusters[ic]) continue;
 594         Float_t y = yexp[iTime] - clusters[ic]->GetY();
 595         Float_t z = zexp[iTime] - clusters[ic]->GetZ();
 596         Float_t d = y * y + z * z;
 597         if(d > dold) continue;
 598         dold = d;
 599         iptr = ic;
 600       }
 601     }
 602     fIndexes[iTime]  = chamber->GetTB(iTime)->GetGlobalIndex(iptr);
 603     fClusters[iTime] = clusters[iptr];
 604     fY[iTime]        = clusters[iptr]->GetY();
 605     fZ[iTime]        = clusters[iptr]->GetZ();
 606     lastCluster      = tboundary[iTime];
 607     fN2++;
 608   }
 609
 610   // number of minimum numbers of clusters expected for the tracklet
 611   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
 612   if (fN2 < kClmin){
 613     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
 614     fN2 = 0;
 615     return kFALSE;
 616   }
 617
 618   // update used clusters
 619   fNUsed = 0;
 620   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 621     if(!fClusters[iTime]) continue;
 622     if((fClusters[iTime]->IsUsed())) fNUsed++;
 623   }
 624
 625   if (fN2-fNUsed < kClmin){
 626     AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
 627     fN2 = 0;
 628     return kFALSE;
 629   }
 630
 631   return kTRUE;
 632 }
 633
 634 //____________________________________________________________
 635 void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
 636 {
 637 //   Fill in all derived information. It has to be called after recovery from file or HLT.
 638 //   The primitive data are
 639 //   - list of clusters
 640 //   - detector (as the detector will be removed from clusters)
 641 //   - position of anode wire (fX0) - temporary
 642 //   - track reference position and direction
 643 //   - momentum of the track
 644 //   - time bin length [cm]
 645 //
 646 //   A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
 647 //
 648   fReconstructor = rec;
 649   AliTRDgeometry g;
 650   AliTRDpadPlane *pp = g.GetPadPlane(fDet);
 651   fTilt      = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
 652   fPadLength = pp->GetLengthIPad();
 653   fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
 654   fTgl = fZref[1];
 655   fN = 0; fN2 = 0; fMPads = 0.;
 656   AliTRDcluster **cit = &fClusters[0];
 657   for(Int_t ic = knTimebins; ic--; cit++){
 658     if(!(*cit)) return;
 659     fN++; fN2++;
 660     fX[ic] = (*cit)->GetX() - fX0;
 661     fY[ic] = (*cit)->GetY();
 662     fZ[ic] = (*cit)->GetZ();
 663   }
 664   Update(); // Fit();
 665   CookLabels();
 666   GetProbability();
 667 }
 668
 669
 670 //____________________________________________________________________
 671 Bool_t AliTRDseedV1::Fit(Bool_t tilt)
 672 {
 673   //
 674   // Linear fit of the tracklet
 675   //
 676   // Parameters :
 677   //
 678   // Output :
 679   //  True if successful
 680   //
 681   // Detailed description
 682   // 2. Check if tracklet crosses pad row boundary
 683   // 1. Calculate residuals in the y (r-phi) direction
 684   // 3. Do a Least Square Fit to the data
 685   //
 686
 687   const Int_t kClmin = 8;
 688   const Float_t q0 = 100.;
 689   const Float_t clSigma0 = 2.E-2;    //[cm]
 690   const Float_t clSlopeQ = -1.19E-2; //[1/cm]
 691
 692   // get track direction
 693   Double_t y0   = fYref[0];
 694   Double_t dydx = fYref[1];
 695   Double_t z0   = fZref[0];
 696   Double_t dzdx = fZref[1];
 697   Double_t yt, zt;
 698
 699   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
 700   AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
 701
 702   // convertion factor from square to gauss distribution for sigma
 703   Double_t convert = 1./TMath::Sqrt(12.);
 704
 705   // book cluster information
 706   Double_t xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins], sz[knTimebins];
 707   Int_t zRow[knTimebins];
 708
 709
 710   fN = 0;
 711   AliTRDcluster *c=0x0, **jc = &fClusters[0];
 712   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
 713     zRow[ic] = -1;
 714     xc[ic]  = -1.;
 715     yc[ic]  = 999.;
 716     zc[ic]  = 999.;
 717     sy[ic]  = 0.;
 718     sz[ic]  = 0.;
 719     if(!(c = (*jc))) continue;
 720     if(!c->IsInChamber()) continue;
 721     Float_t w = 1.;
 722     if(c->GetNPads()>4) w = .5;
 723     if(c->GetNPads()>5) w = .2;
 724     zRow[fN] = c->GetPadRow();
 725     xc[fN]   = fX0 - c->GetX();
 726     yc[fN]   = c->GetY();
 727     zc[fN]   = c->GetZ();
 728
 729     // extrapolated y value for the track
 730     yt = y0 - xc[fN]*dydx;
 731     // extrapolated z value for the track
 732     zt = z0 - xc[fN]*dzdx;
 733     // tilt correction
 734     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt);
 735
 736     // elaborate cluster error
 737     Float_t qr = c->GetQ() - q0;
 738     sy[fN]   = qr < 0. ? clSigma0*TMath::Exp(clSlopeQ*qr) : clSigma0;
 739
 740     fitterY.AddPoint(&xc[fN], yc[fN]-yt, sy[fN]);
 741
 742     sz[fN]   = fPadLength*convert;
 743     fitterZ.AddPoint(&xc[fN], zc[fN], sz[fN]);
 744     fN++;
 745   }
 746   // to few clusters
 747   if (fN < kClmin) return kFALSE;
 748
 749   // fit XY plane
 750   fitterY.Eval();
 751   fYfit[0] = y0+fitterY.GetFunctionParameter(0);
 752   fYfit[1] = dydx-fitterY.GetFunctionParameter(1);
 753
 754   // check par row crossing
 755   Int_t zN[2*AliTRDseed::knTimebins];
 756   Int_t nz = AliTRDtrackerV1::Freq(fN, zRow, zN, kFALSE);
 757   // more than one pad row crossing
 758   if(nz>2) return kFALSE;
 759
 760
 761   // determine z offset of the fit
 762   Float_t zslope = 0.;
 763   Int_t nchanges = 0, nCross = 0;
 764   if(nz==2){ // tracklet is crossing pad row
 765     // Find the break time allowing one chage on pad-rows
 766     // with maximal number of accepted clusters
 767     Int_t padRef = zRow[0];
 768     for (Int_t ic=1; ic<fN; ic++) {
 769       if(zRow[ic] == padRef) continue;
 770
 771       // debug
 772       if(zRow[ic-1] == zRow[ic]){
 773         printf("ERROR in pad row change!!!\n");
 774       }
 775
 776       // evaluate parameters of the crossing point
 777       Float_t sx = (xc[ic-1] - xc[ic])*convert;
 778       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
 779       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
 780       fCross[3] = TMath::Max(dzdx * sx, .01);
 781       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
 782       padRef    = zRow[ic];
 783       nCross    = ic;
 784       nchanges++;
 785     }
 786   }
 787
 788   // condition on nCross and reset nchanges TODO
 789
 790   if(nchanges==1){
 791     if(dzdx * zslope < 0.){
 792       AliInfo("tracklet direction does not correspond to the track direction. TODO.");
 793     }
 794     SetBit(kRowCross, kTRUE); // mark pad row crossing
 795     fitterZ.AddPoint(&fCross[0], fCross[2], fCross[3]);
 796     fitterZ.Eval();
 797     //zc[nc] = fitterZ.GetFunctionParameter(0);
 798     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
 799     fCross[0] = fX0 - fCross[0];
 800   } else if(nchanges > 1){ // debug
 801     AliError("N pad row crossing > 1.");
 802     return kFALSE;
 803   }
 804
 805   UpdateUsed();
 806
 807   return kTRUE;
 808 }
 809
 810
 811 //___________________________________________________________________
 812 void AliTRDseedV1::Print(Option_t *o) const
 813 {
 814   //
 815   // Printing the seedstatus
 816   //
 817
 818   AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
 819   AliInfo(Form("Nattach[%2d] Nfit[%2d] Nuse[%2d] pads[%f]", fN, fN2, fNUsed, fMPads));
 820   AliInfo(Form("x[%7.2f] y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fX0, fYfit[0], fZfit[0], fYfit[1], fZfit[1]));
 821   AliInfo(Form("Ref        y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
 822
 823
 824   if(strcmp(o, "a")!=0) return;
 825
 826   AliTRDcluster* const* jc = &fClusters[0];
 827   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
 828     if(!(*jc)) continue;
 829     (*jc)->Print(o);
 830   }
 831
 832 /*  printf("  fSigmaY =%f\n", fSigmaY);
 833   printf("  fSigmaY2=%f\n", fSigmaY2);
 834   printf("  fMeanz  =%f\n", fMeanz);
 835   printf("  fZProb  =%f\n", fZProb);
 836   printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);*/
 837
 838 /*  printf("  fC      =%f\n", fC);
 839   printf("  fCC     =%f\n",fCC);
 840   printf("  fChi2   =%f\n", fChi2);
 841   printf("  fChi2Z  =%f\n", fChi2Z);*/
 842 }
 843
 844
 845 //___________________________________________________________________
 846 Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
 847 {
 848   // Checks if current instance of the class has the same essential members
 849   // as the given one
 850
 851   if(!o) return kFALSE;
 852   const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
 853   if(!inTracklet) return kFALSE;
 854
 855   for (Int_t i = 0; i < 2; i++){
 856     if ( fYref[i] != inTracklet->GetYref(i) ) return kFALSE;
 857     if ( fZref[i] != inTracklet->GetZref(i) ) return kFALSE;
 858   }
 859
 860   if ( fSigmaY != inTracklet->GetSigmaY() ) return kFALSE;
 861   if ( fSigmaY2 != inTracklet->GetSigmaY2() ) return kFALSE;
 862   if ( fTilt != inTracklet->GetTilt() ) return kFALSE;
 863   if ( fPadLength != inTracklet->GetPadLength() ) return kFALSE;
 864
 865   for (Int_t i = 0; i < knTimebins; i++){
 866     if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
 867     if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
 868     if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
 869     if ( fIndexes[i] != inTracklet->GetIndexes(i) ) return kFALSE;
 870     if ( fUsable[i] != inTracklet->IsUsable(i) ) return kFALSE;
 871   }
 872
 873   for (Int_t i=0; i < 2; i++){
 874     if ( fYfit[i] != inTracklet->GetYfit(i) ) return kFALSE;
 875     if ( fZfit[i] != inTracklet->GetZfit(i) ) return kFALSE;
 876     if ( fYfitR[i] != inTracklet->GetYfitR(i) ) return kFALSE;
 877     if ( fZfitR[i] != inTracklet->GetZfitR(i) ) return kFALSE;
 878     if ( fLabels[i] != inTracklet->GetLabels(i) ) return kFALSE;
 879   }
 880
 881   if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
 882   if ( fZProb != inTracklet->GetZProb() ) return kFALSE;
 883   if ( fN2 != inTracklet->GetN2() ) return kFALSE;
 884   if ( fNUsed != inTracklet->GetNUsed() ) return kFALSE;
 885   if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
 886   if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
 887   if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
 888
 889   if ( fC != inTracklet->GetC() ) return kFALSE;
 890   if ( fCC != inTracklet->GetCC() ) return kFALSE;
 891   if ( fChi2 != inTracklet->GetChi2() ) return kFALSE;
 892   //  if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
 893
 894   if ( fDet != inTracklet->GetDetector() ) return kFALSE;
 895   if ( fMom != inTracklet->GetMomentum() ) return kFALSE;
 896   if ( fdX != inTracklet->GetdX() ) return kFALSE;
 897
 898   for (Int_t iCluster = 0; iCluster < knTimebins; iCluster++){
 899     AliTRDcluster *curCluster = fClusters[iCluster];
 900     AliTRDcluster *inCluster = inTracklet->GetClusters(iCluster);
 901     if (curCluster && inCluster){
 902       if (! curCluster->IsEqual(inCluster) ) {
 903         curCluster->Print();
 904         inCluster->Print();
 905         return kFALSE;
 906       }
 907     } else {
 908       // if one cluster exists, and corresponding
 909       // in other tracklet doesn't - return kFALSE
 910       if(curCluster || inCluster) return kFALSE;
 911     }
 912   }
 913   return kTRUE;
 914 }