TRD/AliTRDtrackerFitter.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 //  Track fitter                                                             //
  21 //                                                                           //
  22 //  Authors:                                                                 //
  23 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
  24 //    Markus Fasel <M.Fasel@gsi.de>                                          //
  25 //                                                                           //
  26 ///////////////////////////////////////////////////////////////////////////////
  27
  28 #include <TROOT.h>
  29 #include <TLinearFitter.h>
  30 #include <TMath.h>
  31 #include <TTreeStream.h>
  32
  33 #include "AliLog.h"
  34
  35 #include "AliTRDseed.h"
  36 #include "AliTRDcalibDB.h"
  37 #include "AliTRDcluster.h"
  38 #include "AliTRDReconstructor.h"
  39 #include "AliTRDseedV1.h"
  40 #include "AliTRDtrackerFitter.h"
  41
  42 #define DEBUG
  43
  44 ClassImp(AliTRDtrackerFitter)
  45
  46 //_________________________________________________________________________
  47 AliTRDtrackerFitter::AliTRDtrackerFitter()
  48   :TObject()
  49   ,fFitterTC(0x0)
  50   ,fFitterT2(0x0)
  51   ,fRieman1(0x0)
  52   ,fRieman2(0x0)
  53   ,fChi2TR(0)
  54   ,fChi2TC(0)
  55   ,fCR(0)
  56   ,fCC(0)
  57   ,fDca(0)
  58   ,fDzmf(0)
  59   ,fZmf(0)
  60   ,fNlayers(0)
  61   ,fDebugStream(0x0)
  62 {
  63   //
  64   // Constructor
  65   //
  66
  67         fRieman1   = new AliRieman(1000);
  68         fRieman2  = new AliRieman(1000);
  69         fFitterTC = new TLinearFitter(2,"hyp2");
  70         fFitterT2 = new TLinearFitter(4,"hyp4");
  71         fFitterTC->StoreData(kTRUE);
  72         fFitterT2->StoreData(kTRUE);
  73 }
  74
  75 //_________________________________________________________________________
  76 AliTRDtrackerFitter::AliTRDtrackerFitter(const AliTRDtrackerFitter &f)
  77   :TObject(f)
  78   ,fFitterTC(0x0)
  79   ,fFitterT2(0x0)
  80   ,fRieman1(0x0)
  81   ,fRieman2(0x0)
  82   ,fChi2TR(f.fChi2TR)
  83   ,fChi2TC(f.fChi2TC)
  84   ,fCR(f.fCR)
  85   ,fCC(f.fCC)
  86   ,fDca(f.fDca)
  87   ,fDzmf(f.fDzmf)
  88   ,fZmf(f.fZmf)
  89   ,fNlayers(f.fNlayers)
  90   ,fDebugStream(0x0)
  91 {
  92   //
  93   // Copy Constructor (performs a deep copy)
  94   //
  95
  96         fRieman1  = new AliRieman(*f.fRieman1);
  97         fRieman2  = new AliRieman(*f.fRieman2);
  98         fFitterTC = new TLinearFitter(*f.fFitterTC);
  99         fFitterT2 = new TLinearFitter(*f.fFitterT2);
 100         fFitterTC->StoreData(kTRUE);
 101         fFitterT2->StoreData(kTRUE);
 102 }
 103
 104 //_________________________________________________________________________
 105 AliTRDtrackerFitter::~AliTRDtrackerFitter()
 106 {
 107   //
 108   // Destructor
 109   //
 110
 111         delete fRieman1;
 112         delete fRieman2;
 113         delete fFitterTC;
 114         delete fFitterT2;
 115 }
 116
 117 //_________________________________________________________________________
 118 AliTRDtrackerFitter &AliTRDtrackerFitter::operator=(const AliTRDtrackerFitter& fitter)
 119 {
 120   //
 121   // Assignment operator using the copy Method of this class
 122   //
 123         if(this != &fitter){
 124                 fitter.Copy(*this);
 125         }
 126         return *this;
 127 }
 128
 129 //_________________________________________________________________________
 130 void AliTRDtrackerFitter::Copy(TObject &f) const
 131 {
 132   //
 133   // Copy method.
 134   // Performs a deep copy. Mainly used in the Assignment operator
 135   //
 136
 137   AliTRDtrackerFitter &tf = (AliTRDtrackerFitter &) f;
 138
 139         if(tf.fFitterTC)
 140                 delete tf.fFitterTC;
 141         tf.fFitterTC = new TLinearFitter(*fFitterTC);
 142         if(tf.fFitterT2)
 143                 delete tf.fFitterT2;
 144         tf.fFitterT2 = new TLinearFitter(*fFitterT2);
 145         if(tf.fRieman1)
 146                 delete tf.fRieman1;
 147         tf.fRieman1 = new AliRieman(*fRieman1);
 148         if(tf.fRieman2)
 149                 delete fRieman2;
 150         tf.fRieman2 = new AliRieman(*fRieman2);
 151         tf.fChi2TR = fChi2TR;
 152         tf.fChi2TC = fChi2TC;
 153         tf.fCR = fCR;
 154         tf.fCC = fCC;
 155         tf.fDca =fDca;
 156         tf.fDzmf = fDzmf;
 157         tf.fZmf =fZmf;
 158         tf.fNlayers = fNlayers;
 159         tf.fDebugStream = 0x0;
 160         fFitterTC->StoreData(kTRUE);
 161         fFitterT2->StoreData(kTRUE);
 162 }
 163
 164 //_________________________________________________________________________
 165 void AliTRDtrackerFitter::FitRieman(AliTRDcluster **cl, Int_t nLayers)
 166 {
 167   //
 168   // Performs a Rieman Fit including 4 clusters (one in each layer)
 169   //
 170         fRieman1->Reset();
 171         for(Int_t i = 0; i < nLayers; i++)
 172                 fRieman1->AddPoint(cl[i]->GetX(), cl[i]->GetY(), cl[i]->GetZ(), 1, 10);
 173         fRieman1->Update();
 174 }
 175
 176 //_________________________________________________________________________
 177 void AliTRDtrackerFitter::FitRieman(AliTRDseedV1 *cseed, Int_t *planes)
 178 {
 179   //
 180   // Performs a Rieman fit using four respectively six seeds
 181   // 2 times used: 1st with 4 parameters and Later in the full track
 182   // fitting Part with 6 parameters
 183   //
 184
 185         Int_t lplanes[] = {0, 1, 2, 3, 4, 5}, *tplanes;
 186         Int_t nplanes;
 187         if(planes){
 188                 nplanes = 4;
 189                 tplanes = planes;
 190         } else {
 191                 nplanes = 6;
 192                 tplanes = &lplanes[0];
 193         }
 194         fRieman1->Reset();
 195         for(Int_t iLayer = 0; iLayer < nplanes; iLayer++){
 196                 if(!cseed[tplanes[iLayer]].IsOK()) continue;
 197                 fRieman1->AddPoint(cseed[tplanes[iLayer]].GetX0(), cseed[tplanes[iLayer]].GetYfitR(0), cseed[tplanes[iLayer]].GetZProb(), 1, 10);
 198         }
 199         fRieman1->Update();
 200 }
 201
 202 //_________________________________________________________________________
 203 Double_t AliTRDtrackerFitter::FitHyperplane(AliTRDseedV1 *cseed
 204                                           , Double_t chi2ZF
 205                                           , Double_t zFitter)
 206 {
 207   //
 208   // performs a hyperplane fit
 209   //
 210   // Two linear fitters:  one in which kz is fixed and one in which kz is a free parameter
 211   // checking if values are acceptable and improves the fitresults
 212   // Calculates curvature parameters and chisquares
 213   // Return: likelihood value as a measurement for the seedquality
 214   //
 215
 216   //const Int_t kMaxTimebins = 100; // Buffer
 217         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 218         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 219         Int_t npointsT = 0;
 220         fFitterTC->ClearPoints();
 221         fFitterT2->ClearPoints();
 222
 223         Double_t fHl[6];
 224         Double_t xref2 = (cseed[2].GetX0() + cseed[3].GetX0())/2;
 225         fRieman2->Reset();
 226         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 227                 if (!cseed[iLayer].IsOK()) continue;
 228                 fHl[iLayer] = cseed[iLayer].GetTilt();
 229
 230                 for (Int_t itime = 0; itime < nTimeBins; itime++) {
 231                         if (!cseed[iLayer].IsUsable(itime)) continue;
 232
 233                         // X relative to the middle chamber
 234                         Double_t x  = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
 235                         Double_t y  = cseed[iLayer].GetY(itime);
 236                         Double_t z  = cseed[iLayer].GetZ(itime);
 237                         // ExB correction to the correction
 238                         // Tilted fRieman1
 239                         Double_t uvt[6];
 240                         // Global x
 241                         Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
 242                         Double_t t  = 1.0 / (x2*x2 + y*y);
 243                         uvt[1] = t;                 // t
 244                         uvt[0] = 2.0 * x2 * uvt[1]; // u
 245                         uvt[2] = 2.0 * fHl[iLayer] * uvt[1];
 246                         uvt[3] = 2.0 * fHl[iLayer] * x * uvt[1];
 247                         uvt[4] = 2.0 * (y + fHl[iLayer]*z) * uvt[1];
 248                         Double_t error = 2.0 * 0.2 * uvt[1];
 249                         fFitterT2->AddPoint(uvt,uvt[4],error);
 250
 251                         //Constrained fRieman1
 252                         z = cseed[iLayer].GetZ(itime);
 253                         uvt[0] = 2.0 * x2 * t; // u
 254                         uvt[1] = 2.0 * fHl[iLayer] * x2 * uvt[1];
 255                         uvt[2] = 2.0 * (y + fHl[iLayer] * (z - zFitter)) * t;   // parameter that is most propably wrong
 256                         fFitterTC->AddPoint(uvt,uvt[2],error);
 257                         fRieman2->AddPoint(x2,y,z,1,10);
 258                         npointsT++;
 259                 }
 260         } // Loop: iLayer
 261
 262         fRieman2->Update();
 263         fFitterTC->Eval();
 264         fFitterT2->Eval();
 265         Double_t rpolz0 = fFitterT2->GetParameter(3);
 266         Double_t rpolz1 = fFitterT2->GetParameter(4);
 267
 268         // Linear fitter  - not possible to make boundaries
 269         // Do not accept non possible z and dzdx combinations
 270         Bool_t acceptablez = kTRUE;
 271         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 272                 if(!cseed[iLayer].IsOK()) continue;
 273                 Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
 274                 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > cseed[iLayer].GetPadLength() * 0.5 + 1.0) acceptablez = kFALSE;
 275         }
 276         if (!acceptablez) {
 277                 fFitterT2->FixParameter(3,fRieman1->GetZat(xref2));
 278                 fFitterT2->FixParameter(4,fRieman1->GetDZat(xref2));
 279                 fFitterT2->Eval();
 280                 fFitterT2->ReleaseParameter(3);
 281                 fFitterT2->ReleaseParameter(4);
 282                 rpolz0 = fFitterT2->GetParameter(3);
 283                 rpolz1 = fFitterT2->GetParameter(4);
 284         }
 285
 286         fChi2TR = fFitterT2->GetChisquare() / Float_t(npointsT);
 287         fChi2TC = fFitterTC->GetChisquare() / Float_t(npointsT);
 288         Double_t polz1c = fFitterTC->GetParameter(2);
 289         Double_t polz0c = polz1c * xref2;
 290         Double_t aC     =  fFitterTC->GetParameter(0);
 291         Double_t bC     =  fFitterTC->GetParameter(1);
 292         fCC     =  aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
 293         Double_t aR     =  fFitterT2->GetParameter(0);
 294         Double_t bR     =  fFitterT2->GetParameter(1);
 295         Double_t dR     =  fFitterT2->GetParameter(2);
 296         fCR     =  1.0 + bR*bR - dR*aR;
 297         fDca    =  0.0;
 298         if (fCR > 0.0) {
 299                 fDca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
 300                 fCR  =  aR / TMath::Sqrt(fCR);
 301         }
 302
 303
 304         Double_t chi2ZT2 = 0.0;
 305         Double_t chi2ZTC = 0.0;
 306         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 307                 if(!cseed[iLayer].IsOK()) continue;
 308                 Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
 309                 Double_t zTC = polz0c + polz1c * (cseed[iLayer].GetX0() - xref2);
 310                 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
 311                 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
 312         }
 313         chi2ZT2 /= TMath::Max((fNlayers - 3.0),1.0);
 314         chi2ZTC /= TMath::Max((fNlayers - 3.0),1.0);
 315
 316
 317         AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
 318         Float_t sumdaf = 0.0;
 319         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 320                 if(!cseed[iLayer].IsOK()) continue;
 321                 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))/ cseed[iLayer].GetSigmaY2());
 322         }
 323         sumdaf /= Float_t (fNlayers - 2.0);
 324
 325         // Cook Likelihoods
 326         Double_t likezf     = TMath::Exp(-chi2ZF * 0.14);
 327         Double_t likechi2C  = TMath::Exp(-fChi2TC * 0.677);
 328         Double_t likechi2TR = TMath::Exp(-fChi2TR * 0.78);
 329         Double_t likeaf     = TMath::Exp(-sumdaf * 3.23);
 330         Double_t likef = likezf * likechi2TR * likeaf;
 331
 332 #ifdef DEBUG
 333         if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 2){
 334                 TTreeSRedirector &treeStreamer = *fDebugStream;
 335                 treeStreamer << "FitHyperplane"
 336                         << "seed0.="     << &cseed[0]
 337                         << "seed1.="     << &cseed[1]
 338                         << "seed2.="     << &cseed[2]
 339                         << "seed3.="     << &cseed[3]
 340                         << "seed4.="     << &cseed[4]
 341                         << "seed5.="     << &cseed[5]
 342                         << "chi2TR="     << fChi2TR
 343                         << "chi2TC="     << fChi2TC
 344                         << "chi2ZT2="    << chi2ZT2
 345                         << "chi2ZTC="    << chi2ZTC
 346                         << "CC="         << fCC
 347                         << "CR="         << fCR
 348                         << "DR="         << dR
 349                         << "DCA="        << fDca
 350                         << "Polz0="      << polz0c
 351                         << "Polz1="      << polz1c
 352                         << "RPolz0="     << rpolz0
 353                         << "RPolz1="     << rpolz1
 354                         << "Likechi2C="  << likechi2C
 355                         << "Likechi2TR=" << likechi2TR
 356                         << "Likezf="     << likezf
 357                         << "Likeaf="     << likeaf
 358                         << "LikeF="      << likef
 359                         << "Rieman1.="   << fRieman1
 360                         << "Rieman2.="   << fRieman2
 361                         << "\n";
 362         }
 363 #endif
 364
 365         return likef;
 366 }
 367
 368 //_________________________________________________________________________
 369 void AliTRDtrackerFitter::GetHyperplaneFitChi2(Double_t *chisquares) const
 370 {
 371   //
 372   // Getter method returns the chisquares of the hyperplane fit
 373   //
 374
 375         chisquares[0] = fChi2TR;
 376         chisquares[1] = fChi2TC;
 377 }
 378
 379 //_________________________________________________________________________
 380 void AliTRDtrackerFitter::GetHyperplaneFitResults(Double_t *params) const
 381 {
 382   //
 383   // Getter Method returning the Curvature parameters of the hyperplane fit
 384   //
 385
 386         params[0] = fCC;
 387         params[1] = fCR;
 388         params[2] = fDca;
 389 }
 390
 391 //_________________________________________________________________________
 392 void AliTRDtrackerFitter::Reset()
 393 {
 394   //
 395   // Resets the object
 396   //
 397
 398         fRieman1->Reset();
 399         fRieman2->Reset();
 400         //fFitterTC->Clear();
 401         //fFitterT2->Clear();
 402 }