TPC/TPCcalib/AliTPCcalibTracks.cxx

   1
   2 #include <TROOT.h>
   3 #include <TChain.h>
   4 #include <TFile.h>
   5 #include <TH3F.h>
   6 //
   7 #include <TPDGCode.h>
   8 #include <TStyle.h>
   9 #include "TLinearFitter.h"
  10 #include "TMatrixD.h"
  11 #include "TTreeStream.h"
  12 #include "TF1.h"
  13
  14
  15
  16 #include "AliMagF.h"
  17 #include "AliTracker.h"
  18 #include "AliESD.h"
  19 #include "AliESDtrack.h"
  20 #include "AliESDfriend.h"
  21 #include "AliESDfriendTrack.h"
  22 #include "AliTPCseed.h"
  23 #include "AliTPCclusterMI.h"
  24 #include "AliTPCROC.h"
  25
  26
  27 #include "AliTPCParamSR.h"
  28 #include "AliTPCClusterParam.h"
  29 #include "AliTrackPointArray.h"
  30 #include "TCint.h"
  31 #include "AliTPCcalibTracks.h"
  32
  33 ClassImp(AliTPCcalibTracks)
  34
  35 AliTPCParam  param;
  36
  37
  38
  39 AliTPCcalibTracks::AliTPCcalibTracks() :
  40    TNamed(),
  41    fHclus(0),
  42    fFileNo(0)
  43  {
  44    G__SetCatchException(0);
  45    param.Update();
  46    TFile fparam("/u/miranov/TPCClusterParam.root");
  47    fClusterParam  =  (AliTPCClusterParam *) fparam.Get("Param");
  48    if (fClusterParam){
  49      //fClusterParam->SetInstance(fClusterParam);
  50    }else{
  51      printf("Cluster Param not found\n");
  52    }
  53    fDebugStream = new TTreeSRedirector("TPCSelectorDebug.root");
  54  }
  55
  56
  57 void    AliTPCcalibTracks::ProcessTrack(AliTPCseed * seed){
  58 }
  59
  60
  61 Int_t   AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
  62   //
  63   // calculate bins for given q and pad type
  64   // used in TObjArray
  65   //
  66   Int_t   res  = TMath::Max(TMath::Nint((TMath::Sqrt(q)-3.)),0);
  67   res*=3;
  68   res+=pad;
  69   return res;
  70 }
  71
  72 Int_t   AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
  73   //
  74   // calculate bins for given iq and pad type
  75   // used in TObjArray
  76   //
  77   return iq*3+pad;;
  78 }
  79
  80 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
  81   Int_t bin0 = bin/3;
  82   bin0+=3;
  83   return bin0*bin0;
  84 }
  85
  86
  87
  88
  89
  90
  91
  92 void AliTPCcalibTracks::ProofSlaveBegin(TList * output)
  93 {
  94   // Called on PROOF - fill output list
  95   //fChain = tree;
  96   //Init(tree);
  97
  98   char chname[1000];
  99   TProfile * prof1=0;
 100   TH1F     * his1 =0;
 101   fHclus = new TH1I("hclus","Number of clusters",100,0,200);
 102   output->AddLast(fHclus);
 103
 104
 105   //
 106   // Amplitude  - sector -row histograms
 107   //
 108   fArrayAmpRow = new TObjArray(72);
 109   fArrayAmp    = new TObjArray(72);
 110   for (Int_t i=0; i<36; i++){
 111     sprintf(chname,"Amp_row_Sector%d",i);
 112     prof1 = new TProfile(chname,chname,63,0,64);
 113     prof1->SetXTitle("Pad row");
 114     prof1->SetYTitle("Mean Max amplitude");
 115     fArrayAmpRow->AddAt(prof1,i);
 116     output->AddLast(prof1);
 117     sprintf(chname,"Amp_row_Sector%d",i+36);
 118     prof1 = new TProfile(chname,chname,96,0,97);
 119     prof1->SetXTitle("Pad row");
 120     prof1->SetYTitle("Mean Max  amplitude");
 121     fArrayAmpRow->AddAt(prof1,i+36);
 122     output->AddLast(prof1);
 123     //
 124     // amplitude
 125     sprintf(chname,"Amp_Sector%d",i);
 126     his1 = new TH1F(chname,chname,250,0,500);
 127     his1->SetXTitle("Max Amplitude (ADC)");
 128     fArrayAmp->AddAt(his1,i);
 129     output->AddLast(his1);
 130     sprintf(chname,"Amp_Sector%d",i+36);
 131     his1 = new TH1F(chname,chname,200,0,600);
 132     his1->SetXTitle("Max Amplitude (ADC)");
 133     fArrayAmp->AddAt(his1,i+36);
 134     output->AddLast(his1);
 135     //
 136   }
 137
 138   fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
 139   fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
 140   output->AddLast(fDeltaY);
 141   output->AddLast(fDeltaZ);
 142
 143   fResolY = new TObjArray(3);
 144   fResolZ = new TObjArray(3);
 145   fRMSY   = new TObjArray(3);
 146   fRMSZ   = new TObjArray(3);
 147   TH3F * his3D;
 148   //
 149   his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
 150   fResolY->AddAt(his3D,0);
 151   output->AddLast(his3D);
 152   his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
 153   fResolY->AddAt(his3D,1);
 154   output->AddLast(his3D);
 155   his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
 156   fResolY->AddAt(his3D,2);
 157   output->AddLast(his3D);
 158   //
 159   his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
 160   fResolZ->AddAt(his3D,0);
 161   output->AddLast(his3D);
 162   his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
 163   fResolZ->AddAt(his3D,1);
 164   output->AddLast(his3D);
 165   his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
 166   fResolZ->AddAt(his3D,2);
 167   output->AddLast(his3D);
 168   //
 169   his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
 170   fRMSY->AddAt(his3D,0);
 171   output->AddLast(his3D);
 172   his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
 173   fRMSY->AddAt(his3D,1);
 174   output->AddLast(his3D);
 175   his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
 176   fRMSY->AddAt(his3D,2);
 177   output->AddLast(his3D);
 178   //
 179   his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
 180   fRMSZ->AddAt(his3D,0);
 181   output->AddLast(his3D);
 182   his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
 183   fRMSZ->AddAt(his3D,1);
 184   output->AddLast(his3D);
 185   his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
 186   fRMSZ->AddAt(his3D,2);
 187   output->AddLast(his3D);
 188   //
 189   fArrayQDY = new TObjArray(300);
 190   fArrayQDZ = new TObjArray(300);
 191   fArrayQRMSY = new TObjArray(300);
 192   fArrayQRMSZ = new TObjArray(300);
 193   for (Int_t iq=0; iq<10; iq++){
 194     for (Int_t ipad=0; ipad<3; ipad++){
 195       Int_t   bin   = GetBin(iq,ipad);
 196       Float_t qmean = GetQ(bin);
 197       char name[200];
 198       sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean);
 199       his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
 200       fArrayQDY->AddAt(his3D,bin);
 201       output->AddLast(his3D);
 202       sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
 203       his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
 204       fArrayQDZ->AddAt(his3D,bin);
 205       output->AddLast(his3D);
 206       //
 207       sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean);
 208       his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
 209       fArrayQRMSY->AddAt(his3D,bin);
 210       output->AddLast(his3D);
 211       sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
 212       his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
 213       fArrayQRMSZ->AddAt(his3D,bin);
 214       output->AddLast(his3D);
 215     }
 216   }
 217 }
 218
 219
 220
 221
 222 Float_t AliTPCcalibTracks::TPCBetheBloch(Float_t bg)
 223 {
 224  //
 225  // Bethe-Bloch energy loss formula
 226  //
 227  const Double_t kp1=0.76176e-1;
 228  const Double_t kp2=10.632;
 229  const Double_t kp3=0.13279e-4;
 230  const Double_t kp4=1.8631;
 231  const Double_t kp5=1.9479;
 232  Double_t dbg = (Double_t) bg;
 233  Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
 234  Double_t aa = TMath::Power(beta,kp4);
 235  Double_t bb = TMath::Power(1./dbg,kp5);
 236  bb=TMath::Log(kp3+bb);
 237  return ((Float_t)((kp2-aa-bb)*kp1/aa));
 238 }
 239
 240
 241 Bool_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
 242   //
 243   //
 244   //
 245   const Int_t   kMinClusters  = 20;
 246   const Float_t kMinRatio     = 0.4;
 247   const Float_t kMax1pt       = 0.5;
 248   const Float_t kEdgeYXCutNoise    = 0.13;
 249   const Float_t kEdgeThetaCutNoise = 0.018;
 250   //
 251   // edge induced noise tracks - NEXT RELEASE will be removed during tracking
 252   if (TMath::Abs(track->GetY()/track->GetX())> kEdgeYXCutNoise)
 253     if (TMath::Abs(track->GetTgl())<kEdgeThetaCutNoise) return kFALSE;
 254
 255   //
 256   if (track->GetNumberOfClusters()<kMinClusters) return kFALSE;
 257   Float_t ratio = track->GetNumberOfClusters()/(track->GetNFoundable()+1.);
 258   if (ratio<kMinRatio) return kFALSE;
 259   Float_t mpt = track->Get1Pt();
 260   if (TMath::Abs(mpt)>kMax1pt) return kFALSE;
 261   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
 262   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
 263   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
 264
 265   return kTRUE;
 266 }
 267
 268 void AliTPCcalibTracks::FillHistoCluster(AliTPCseed * track){
 269   //
 270   //
 271   //
 272   const Int_t kFirstLargePad = 127;
 273   const Float_t kLargePadSize = 1.5;
 274   for (Int_t irow=0; irow<159; irow++){
 275     AliTPCclusterMI * cluster = track->GetClusterPointer(irow);
 276     if (!cluster) continue;
 277     Int_t sector = cluster->GetDetector();
 278     if (cluster->GetQ()<=0) continue;
 279     Float_t max = cluster->GetMax();
 280     printf ("irow, kFirstLargePad = %d, %d \n",irow,kFirstLargePad);
 281     if ( irow >= kFirstLargePad) {
 282       max /= kLargePadSize;
 283     }
 284     TProfile *profAmpRow =  (TProfile*)fArrayAmpRow->At(sector);
 285     profAmpRow->Fill(cluster->GetRow(), max);
 286   }
 287 }
 288
 289 void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
 290   //
 291   // fill resolution histograms - localy - trcklet in the neighborhood
 292   //
 293   const Int_t   kDelta  = 10;      // delta rows to fit
 294   const Float_t kMinRatio = 0.75;  // minimal ratio
 295   const Float_t kCutChi2 = 6.;     // cut chi2 - left right  - kink removal
 296   const Float_t kErrorFraction = 0.5;  // use only clusters with small intrpolation error - for error param
 297   const Int_t   kFirstLargePad = 127;
 298   const Float_t kLargePadSize = 1.5;
 299   static TLinearFitter fitterY2(3,"pol2");
 300   static TLinearFitter fitterZ2(3,"pol2");
 301   static TLinearFitter fitterY0(2,"pol1");
 302   static TLinearFitter fitterZ0(2,"pol1");
 303   static TLinearFitter fitterY1(2,"pol1");
 304   static TLinearFitter fitterZ1(2,"pol1");
 305   TVectorD      paramY0(2);
 306   TVectorD      paramY1(2);
 307   TVectorD      paramY2(3);
 308   TVectorD      paramZ0(2);
 309   TVectorD      paramZ1(2);
 310   TVectorD      paramZ2(3);
 311   TMatrixD    matrixY0(2,2);
 312   TMatrixD    matrixZ0(2,2);
 313   TMatrixD    matrixY1(2,2);
 314   TMatrixD    matrixZ1(2,2);
 315   //
 316   // estimate mean error
 317   //
 318   Int_t nTrackletsAll  = 0;
 319   Int_t nClusters      = 0;
 320   Float_t csigmaY       = 0;
 321   Float_t csigmaZ       = 0;
 322   Int_t sectorG       = -1;
 323   for (Int_t irow=0; irow<159; irow++){
 324     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 325     if (!cluster0) continue;
 326     Int_t sector = cluster0->GetDetector();
 327     if (sector!=sectorG){
 328       nClusters=0;
 329       fitterY2.ClearPoints();
 330       fitterZ2.ClearPoints();
 331       sectorG=sector;
 332     }else{
 333       nClusters++;
 334       Double_t x = cluster0->GetX();
 335       fitterY2.AddPoint(&x,cluster0->GetY(),1);
 336       fitterZ2.AddPoint(&x,cluster0->GetZ(),1);
 337       //
 338       if (nClusters>=kDelta+3){
 339         fitterY2.Eval();
 340         fitterZ2.Eval();
 341         nTrackletsAll++;
 342         csigmaY+=fitterY2.GetChisquare()/(nClusters-3.);
 343         csigmaZ+=fitterZ2.GetChisquare()/(nClusters-3.);
 344         nClusters=-1;
 345         fitterY2.ClearPoints();
 346         fitterZ2.ClearPoints();
 347       }
 348     }
 349   }
 350   csigmaY = TMath::Sqrt(csigmaY/nTrackletsAll);
 351   csigmaZ = TMath::Sqrt(csigmaZ/nTrackletsAll);
 352   //
 353   //
 354   //
 355   for (Int_t irow=0; irow<159; irow++){
 356     Int_t nclFound     = 0;
 357     Int_t nclFoundable = 0;
 358     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 359     if (!cluster0) continue;
 360     Int_t sector = cluster0->GetDetector();
 361     Float_t xref = cluster0->GetX();
 362     //
 363     // check the neighborhood occupancy - (Delta ray - noise removal)
 364     //
 365     for (Int_t idelta= -kDelta; idelta<=kDelta; idelta++){
 366       if (idelta==0) continue;
 367       if (idelta+irow<0) continue;
 368       if (idelta+irow>159) continue;
 369       AliTPCclusterMI * clusterD = track->GetClusterPointer(irow);
 370       if ( clusterD && clusterD->GetDetector()!= sector) continue;
 371       if (clusterD->GetType()<0) continue;
 372       nclFoundable++;
 373       if (clusterD) nclFound++;
 374     }
 375     if (nclFound<kDelta*kMinRatio) continue;
 376     if (Float_t(nclFound)/Float_t(nclFoundable)<kMinRatio) continue;
 377     //
 378     // Make Fit
 379     //
 380     fitterY2.ClearPoints();
 381     fitterZ2.ClearPoints();
 382     fitterY0.ClearPoints();
 383     fitterZ0.ClearPoints();
 384     fitterY1.ClearPoints();
 385     fitterZ1.ClearPoints();
 386     //
 387     nclFound=0;
 388     Int_t ncl0=0;
 389     Int_t ncl1=0;
 390     for (Int_t idelta=-kDelta; idelta<=kDelta; idelta++){
 391       if (idelta==0) continue;
 392       if (idelta+irow<0) continue;
 393       if (idelta+irow>159) continue;
 394       AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
 395       if (!cluster) continue;
 396       if (cluster->GetType()<0) continue;
 397       if (cluster->GetDetector()!=sector) continue;
 398       Double_t x = cluster->GetX()-xref;
 399       nclFound++;
 400       if (idelta<0){
 401         ncl0++;
 402         fitterY0.AddPoint(&x, cluster->GetY(),csigmaY);
 403         fitterZ0.AddPoint(&x, cluster->GetZ(),csigmaZ);
 404       }
 405       if (idelta>0){
 406         ncl1++;
 407         fitterY1.AddPoint(&x, cluster->GetY(),csigmaY);
 408         fitterZ1.AddPoint(&x, cluster->GetZ(),csigmaZ);
 409       }
 410       fitterY2.AddPoint(&x, cluster->GetY(),csigmaY);
 411       fitterZ2.AddPoint(&x, cluster->GetZ(),csigmaZ);
 412     }
 413     if (nclFound<kDelta*kMinRatio) continue;
 414     fitterY2.Eval();
 415     fitterZ2.Eval();
 416     Double_t chi2 = (fitterY2.GetChisquare()+fitterZ2.GetChisquare())/(2.*nclFound-6.);
 417     if (chi2>kCutChi2) continue;
 418     //
 419     //
 420     // REMOVE KINK
 421     //
 422     if (ncl0>4){
 423       fitterY0.Eval();
 424       fitterZ0.Eval();
 425     }
 426     if (ncl1>4){
 427       fitterY1.Eval();
 428       fitterZ1.Eval();
 429     }
 430     //
 431     //
 432     if (ncl0>4&&ncl1>4){
 433       fitterY0.GetCovarianceMatrix(matrixY0);
 434       fitterY1.GetCovarianceMatrix(matrixY1);
 435       fitterZ0.GetCovarianceMatrix(matrixZ0);
 436       fitterZ1.GetCovarianceMatrix(matrixZ1);
 437       fitterY1.GetParameters(paramY1);
 438       fitterZ1.GetParameters(paramZ1);
 439       fitterY0.GetParameters(paramY0);
 440       fitterZ0.GetParameters(paramZ0);
 441       paramY0-= paramY1;
 442       paramZ0-= paramZ1;
 443       matrixY0+= matrixY1;
 444       matrixZ0+= matrixZ1;
 445       Double_t chi2 =0;
 446       TMatrixD difY(2,1,paramY0.GetMatrixArray());
 447       TMatrixD difYT(1,2,paramY0.GetMatrixArray());
 448       matrixY0.Invert();
 449       TMatrixD mulY(matrixY0,TMatrixD::kMult,difY);
 450       TMatrixD chi2Y(difYT,TMatrixD::kMult,mulY);
 451       chi2+=chi2Y(0,0);
 452       TMatrixD difZ(2,1,paramZ0.GetMatrixArray());
 453       TMatrixD difZT(1,2,paramZ0.GetMatrixArray());
 454       matrixZ0.Invert();
 455       TMatrixD mulZ(matrixZ0,TMatrixD::kMult,difZ);
 456       TMatrixD chi2Z(difZT,TMatrixD::kMult,mulZ);
 457       chi2+= chi2Z(0,0);
 458       if (chi2*0.25>kCutChi2) continue;
 459     }
 460     Double_t paramY[4], paramZ[4];
 461     paramY[0] = fitterY2.GetParameter(0);
 462     paramY[1] = fitterY2.GetParameter(1);
 463     paramY[2] = fitterY2.GetParameter(2);
 464     paramZ[0] = fitterZ2.GetParameter(0);
 465     paramZ[1] = fitterZ2.GetParameter(1);
 466     paramZ[2] = fitterZ2.GetParameter(2);
 467     //
 468     //
 469     //
 470     Double_t tracky = paramY[0];
 471     Double_t trackz = paramZ[0];
 472     Float_t  deltay = tracky-cluster0->GetY();
 473     Float_t  deltaz = trackz-cluster0->GetZ();
 474     Float_t  angley = paramY[1]-paramY[0]/xref;
 475     Float_t  anglez = paramZ[1];
 476     //
 477     //
 478     Float_t max = cluster0->GetMax();
 479     TProfile *profAmpRow =  (TProfile*)fArrayAmpRow->At(sector);
 480     if ( irow >= kFirstLargePad) max /= kLargePadSize;
 481     profAmpRow->Fill(cluster0->GetRow(), max);
 482     TH1F *hisAmp =  (TH1F*)fArrayAmp->At(sector);
 483     hisAmp->Fill(max);
 484     //
 485     //
 486     Int_t  ipad=0;
 487     if (cluster0->GetDetector()>=36) {
 488       ipad=1;
 489       if (cluster0->GetRow()>63) ipad=2;
 490     }
 491     //
 492     //
 493     TH3F * his3=0;
 494     his3 = (TH3F*)fRMSY->At(ipad);
 495     if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley),TMath::Sqrt(cluster0->GetSigmaY2()));
 496     his3 = (TH3F*)fRMSZ->At(ipad);
 497     if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez),TMath::Sqrt(cluster0->GetSigmaZ2()));
 498       his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(),ipad));
 499       if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()));
 500       //
 501       his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(),ipad));
 502       if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez),TMath::Sqrt(cluster0->GetSigmaZ2()));
 503
 504     //
 505     // Fill resolution histograms
 506     //
 507     Bool_t useForResol= kTRUE;
 508     if (fitterY2.GetParError(0)>kErrorFraction*csigmaY) useForResol=kFALSE;
 509
 510     if (useForResol){
 511       fDeltaY->Fill(deltay);
 512       fDeltaZ->Fill(deltaz);
 513       his3 = (TH3F*)fResolY->At(ipad);
 514       if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay);
 515       his3 = (TH3F*)fResolZ->At(ipad);
 516       if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz);
 517       his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(),ipad));
 518       if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay);
 519       //
 520       his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(),ipad));
 521       if (his3) his3->Fill(250-TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz);
 522     }
 523     //
 524     //
 525     //
 526     if (useForResol&&nclFound>2*kMinRatio*kDelta){
 527       //
 528       // fill resolution trees
 529       //
 530       static TLinearFitter fitY0(3,"pol2");
 531       static TLinearFitter fitZ0(3,"pol2");
 532       static TLinearFitter fitY2(5,"hyp4");
 533       static TLinearFitter fitZ2(5,"hyp4");
 534       static TLinearFitter fitY2Q(5,"hyp4");
 535       static TLinearFitter fitZ2Q(5,"hyp4");
 536       static TLinearFitter fitY2S(5,"hyp4");
 537       static TLinearFitter fitZ2S(5,"hyp4");
 538       fitY0.ClearPoints();
 539       fitZ0.ClearPoints();
 540       fitY2.ClearPoints();
 541       fitZ2.ClearPoints();
 542       fitY2Q.ClearPoints();
 543       fitZ2Q.ClearPoints();
 544       fitY2S.ClearPoints();
 545       fitZ2S.ClearPoints();
 546
 547       for (Int_t idelta=-kDelta; idelta<=kDelta; idelta++){
 548         if (idelta==0) continue;
 549         if (idelta+irow<0) continue;
 550         if (idelta+irow>159) continue;
 551         AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
 552         if (!cluster) continue;
 553         if (cluster->GetType()<0) continue;
 554         if (cluster->GetDetector()!=sector) continue;
 555         Double_t x = cluster->GetX()-xref;
 556         Double_t sigmaY0 = fClusterParam->GetError0Par(0,ipad,(250.0-TMath::Abs(cluster->GetZ())),TMath::Abs(angley));
 557         Double_t sigmaZ0 = fClusterParam->GetError0Par(1,ipad,(250.0-TMath::Abs(cluster->GetZ())),TMath::Abs(anglez));
 558         //
 559         Double_t sigmaYQ = fClusterParam->GetErrorQPar(0,ipad,(250.0-TMath::Abs(cluster->GetZ())),
 560                                                        TMath::Abs(angley), TMath::Abs(cluster->GetMax()));
 561         Double_t sigmaZQ = fClusterParam->GetErrorQPar(1,ipad,(250.0-TMath::Abs(cluster->GetZ())),
 562                                                        TMath::Abs(anglez),TMath::Abs(cluster->GetMax()));
 563         Double_t sigmaYS = fClusterParam->GetErrorQParScaled(0,ipad,(250.0-TMath::Abs(cluster->GetZ())),
 564                                                        TMath::Abs(angley), TMath::Abs(cluster->GetMax()));
 565         Double_t sigmaZS = fClusterParam->GetErrorQParScaled(1,ipad,(250.0-TMath::Abs(cluster->GetZ())),
 566                                                        TMath::Abs(anglez),TMath::Abs(cluster->GetMax()));
 567         Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,ipad,(250.0-TMath::Abs(cluster->GetZ())),
 568                                                            TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
 569                                                            TMath::Sqrt(cluster0->GetSigmaY2()),0);
 570         Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,ipad,(250.0-TMath::Abs(cluster->GetZ())),
 571                                                            TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
 572                                                            TMath::Sqrt(cluster0->GetSigmaZ2()),0);
 573         sigmaYS  = TMath::Sqrt(sigmaYS*sigmaYS+rmsYFactor*rmsYFactor/12.);
 574         sigmaZS  = TMath::Sqrt(sigmaZS*sigmaZS+rmsZFactor*rmsZFactor/12.+rmsYFactor*rmsYFactor/24.);
 575         //
 576         if (kDelta!=0){
 577           fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
 578           fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
 579         }
 580         Double_t xxx[4];
 581         xxx[0] =  ((idelta+irow)%2==0)? 1:0;
 582         xxx[1] =  x;
 583         xxx[2] =  ((idelta+irow)%2==0)? x:0;
 584         xxx[3] =  x*x;
 585         fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
 586         fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
 587         fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
 588         fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
 589         fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
 590         fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
 591         //
 592       }
 593       //
 594       fitY0.Eval();
 595       fitZ0.Eval();
 596       fitY2.Eval();
 597       fitZ2.Eval();
 598       fitY2Q.Eval();
 599       fitZ2Q.Eval();
 600       fitY2S.Eval();
 601       fitZ2S.Eval();
 602       Float_t chi2Y0 = fitY0.GetChisquare()/(nclFound-3.);
 603       Float_t chi2Z0 = fitZ0.GetChisquare()/(nclFound-3.);
 604       Float_t chi2Y2 = fitY2.GetChisquare()/(nclFound-5.);
 605       Float_t chi2Z2 = fitZ2.GetChisquare()/(nclFound-5.);
 606       Float_t chi2Y2Q = fitY2Q.GetChisquare()/(nclFound-5.);
 607       Float_t chi2Z2Q = fitZ2Q.GetChisquare()/(nclFound-5.);
 608       Float_t chi2Y2S = fitY2S.GetChisquare()/(nclFound-5.);
 609       Float_t chi2Z2S = fitZ2S.GetChisquare()/(nclFound-5.);
 610       //
 611       static  TVectorD    parY0(3);
 612       static  TMatrixD    matY0(3,3);
 613       static  TVectorD    parZ0(3);
 614       static  TMatrixD    matZ0(3,3);
 615       fitY0.GetParameters(parY0);
 616       fitY0.GetCovarianceMatrix(matY0);
 617       fitZ0.GetParameters(parZ0);
 618       fitZ0.GetCovarianceMatrix(matZ0);
 619       //
 620       static  TVectorD    parY2(5);
 621       static  TMatrixD    matY2(5,5);
 622       static  TVectorD    parZ2(5);
 623       static  TMatrixD    matZ2(5,5);
 624       fitY2.GetParameters(parY2);
 625       fitY2.GetCovarianceMatrix(matY2);
 626       fitZ2.GetParameters(parZ2);
 627       fitZ2.GetCovarianceMatrix(matZ2);
 628       //
 629       static  TVectorD    parY2Q(5);
 630       static  TMatrixD    matY2Q(5,5);
 631       static  TVectorD    parZ2Q(5);
 632       static  TMatrixD    matZ2Q(5,5);
 633       fitY2Q.GetParameters(parY2Q);
 634       fitY2Q.GetCovarianceMatrix(matY2Q);
 635       fitZ2Q.GetParameters(parZ2Q);
 636       fitZ2Q.GetCovarianceMatrix(matZ2Q);
 637       static  TVectorD    parY2S(5);
 638       static  TMatrixD    matY2S(5,5);
 639       static  TVectorD    parZ2S(5);
 640       static  TMatrixD    matZ2S(5,5);
 641       fitY2S.GetParameters(parY2S);
 642       fitY2S.GetCovarianceMatrix(matY2S);
 643       fitZ2S.GetParameters(parZ2S);
 644       fitZ2S.GetCovarianceMatrix(matZ2S);
 645       Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
 646       Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
 647       Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
 648       Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
 649       Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
 650       Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
 651       Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
 652       Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
 653       Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
 654       Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
 655       Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
 656       Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
 657       Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
 658       Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
 659       Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
 660       Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
 661       //
 662       // Error parameters
 663       //
 664       Float_t csigmaY0 = fClusterParam->GetError0Par(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
 665       Float_t csigmaZ0 = fClusterParam->GetError0Par(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
 666       //
 667       Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 668                                                      TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
 669       Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 670                                                        TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
 671       Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 672                                                      TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
 673       Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 674                                                        TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
 675       //
 676       // RMS parameters
 677       //
 678       Float_t meanRMSY = 0;
 679       Float_t meanRMSZ = 0;
 680       Int_t   nclRMS=0;
 681       for (Int_t idelta=-2; idelta<=2; idelta++){
 682         if (idelta+irow<0) continue;
 683         if (idelta+irow>159) continue;
 684         AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
 685         if (!cluster) continue;
 686         meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
 687         meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
 688         nclRMS++;
 689       }
 690       meanRMSY /= nclRMS;
 691       meanRMSZ /= nclRMS;
 692
 693       Float_t rmsY     = TMath::Sqrt(cluster0->GetSigmaY2());
 694       Float_t rmsZ     = TMath::Sqrt(cluster0->GetSigmaZ2());
 695       Float_t rmsYT    = fClusterParam->GetRMSQ(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 696                                                 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
 697       Float_t rmsZT    = fClusterParam->GetRMSQ(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 698                                                 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
 699       Float_t rmsYT0    = fClusterParam->GetRMS0(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 700                                                  TMath::Abs(angley));
 701       Float_t rmsZT0    = fClusterParam->GetRMS0(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 702                                                  TMath::Abs(anglez));
 703       Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 704                                                      TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
 705       Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 706                                                      TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
 707       Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 708                                                          TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
 709                                                          rmsY,meanRMSY);
 710       Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 711                                                          TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
 712                                                          rmsZ,meanRMSZ);
 713       //
 714       // cluster debug
 715       //
 716       (*fDebugStream)<<"ResolCl"<<
 717         "Sector="<<sector<<
 718         "Cl.="<<cluster0<<
 719         "CSigmaY0="<<csigmaY0<<   // cluster errorY
 720         "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
 721         "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
 722         "CSigmaZ0="<<csigmaZ0<<   //
 723         "CSigmaZQ="<<csigmaZQ<<
 724         "CSigmaZS="<<csigmaZS<<
 725         "shapeYF="<<rmsYFactor<<
 726         "shapeZF="<<rmsZFactor<<
 727         "rmsY="<<rmsY<<
 728         "rmsZ="<<rmsZ<<
 729         "rmsYM="<<meanRMSY<<
 730         "rmsZM="<<meanRMSZ<<
 731         "rmsYT="<<rmsYT<<
 732         "rmsZT="<<rmsZT<<
 733         "rmsYT0="<<rmsYT0<<
 734         "rmsZT0="<<rmsZT0<<
 735         "rmsYS="<<rmsYSigma<<
 736         "rmsZS="<<rmsZSigma<<
 737         "IPad="<<ipad<<
 738         "Ncl="<<nclFound<<
 739         "PY0.="<<&parY0<<
 740         "PZ0.="<<&parZ0<<
 741         "SigmaY0="<<sigmaY0<<
 742         "SigmaZ0="<<sigmaZ0<<
 743         "\n";
 744       //
 745       // tracklet dubug
 746       //
 747       (*fDebugStream)<<"ResolTr"<<
 748         "IPad="<<ipad<<
 749         "Sector="<<sector<<
 750         "Ncl="<<nclFound<<
 751         "chi2Y0="<<chi2Y0<<
 752         "chi2Z0="<<chi2Z0<<
 753         "chi2Y2="<<chi2Y2<<
 754         "chi2Z2="<<chi2Z2<<
 755         "chi2Y2Q="<<chi2Y2Q<<
 756         "chi2Z2Q="<<chi2Z2Q<<
 757         "chi2Y2S="<<chi2Y2S<<
 758         "chi2Z2S="<<chi2Z2S<<
 759         "PY0.="<<&parY0<<
 760         "PZ0.="<<&parZ0<<
 761         "PY2.="<<&parY2<<
 762         "PZ2.="<<&parZ2<<
 763         "PY2Q.="<<&parY2Q<<
 764         "PZ2Q.="<<&parZ2Q<<
 765         "PY2S.="<<&parY2S<<
 766         "PZ2S.="<<&parZ2S<<
 767         "SigmaY0="<<sigmaY0<<
 768         "SigmaZ0="<<sigmaZ0<<
 769         "SigmaDY0="<<sigmaDY0<<
 770         "SigmaDZ0="<<sigmaDZ0<<
 771         "SigmaY2="<<sigmaY2<<
 772         "SigmaZ2="<<sigmaZ2<<
 773         "SigmaDY2="<<sigmaDY2<<
 774         "SigmaDZ2="<<sigmaDZ2<<
 775         "SigmaY2Q="<<sigmaY2Q<<
 776         "SigmaZ2Q="<<sigmaZ2Q<<
 777         "SigmaDY2Q="<<sigmaDY2Q<<
 778         "SigmaDZ2Q="<<sigmaDZ2Q<<
 779         "SigmaY2S="<<sigmaY2S<<
 780         "SigmaZ2S="<<sigmaZ2S<<
 781         "SigmaDY2S="<<sigmaDY2S<<
 782         "SigmaDZ2S="<<sigmaDZ2S<<
 783         "\n";
 784     }
 785   }
 786 }
 787
 788
 789 void  AliTPCcalibTracks::AlignUpDown(AliTPCseed * track, AliESDtrack * esdTrack){
 790   //
 791   // Make simple parabolic fit
 792   //
 793   const Int_t kMinClusters = 60;
 794   const Int_t kMinClustersSector =15;
 795   const Float_t kSigmaCut = 6;
 796   const Float_t kMaxTan = TMath::Tan(TMath::Pi()*10./180.);
 797   const Float_t kDeadZone = 6.;
 798   const Float_t kMinZ     = 15;
 799   if (track->GetNumberOfClusters()<kMinClusters) return;
 800   if (TMath::Abs(track->GetZ())<kMinZ) return;
 801   //
 802   Int_t nclUp   = 0;
 803   Int_t nclDown = 0;
 804   Int_t rSector =-1;
 805   Float_t refX  = (param.GetInnerRadiusUp()+param.GetOuterRadiusLow())*0.5;
 806   for (Int_t irow=0; irow<159; irow++){
 807     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 808     if (!cluster0) continue;
 809     Int_t sector = cluster0->GetDetector();
 810     if (rSector<0) rSector=sector%36;
 811     if (sector%36 != rSector) continue;
 812     if (  ((TMath::Abs(cluster0->GetY())-kDeadZone)/cluster0->GetX())>kMaxTan) continue;  //remove edge clusters
 813     if (sector>35) nclUp++;
 814     if (sector<36) nclDown++;
 815   }
 816   if (nclUp<kMinClustersSector) return;
 817   if (nclDown<kMinClustersSector) return;
 818   //
 819   //
 820   TLinearFitter fitterY(5,"hyp4");  //fitter with common 2 nd derivation
 821   TLinearFitter fitterZ(5,"hyp4");
 822   //
 823   TLinearFitter fitterY0(3,"pol2");
 824   TLinearFitter fitterZ0(3,"pol2");
 825   TLinearFitter fitterY1(3,"pol2");
 826   TLinearFitter fitterZ1(3,"pol2");
 827   //
 828   Float_t msigmay =1;
 829   Float_t msigmaz =1;
 830   Float_t param0[3];
 831   Float_t param1[3];
 832   Float_t angley=0;
 833   Float_t anglez=0;
 834   //
 835   for (Int_t iter=0; iter<3; iter++){
 836     nclUp  = 0;
 837     nclDown= 0;
 838     for (Int_t irow=0; irow<159; irow++){
 839       AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 840       if (!cluster0) continue;
 841       Int_t sector = cluster0->GetDetector();
 842       if (sector%36 != rSector) continue;
 843       Double_t y = cluster0->GetY();
 844       Double_t z = cluster0->GetZ();
 845       //remove edge clusters
 846       if ( (iter==0) && ((TMath::Abs(cluster0->GetY())-kDeadZone)/cluster0->GetX())>kMaxTan ) continue;
 847       if (iter>0){
 848         Float_t tx = cluster0->GetX()-refX;
 849         Float_t ty = 0;
 850         if (sector<36){
 851           ty = param0[0]+param0[1]*tx+param0[2]*tx*tx;
 852         }else{
 853           ty = param1[0]+param1[1]*tx+param1[2]*tx*tx;
 854         }
 855         if (((TMath::Abs(ty)-kDeadZone)/cluster0->GetX())>kMaxTan) continue;
 856         if (TMath::Abs(ty-y)>kSigmaCut*(msigmay+0.2)) continue;
 857       }
 858       Int_t  ipad=0;
 859       if (cluster0->GetDetector()>=36) {
 860         ipad=1;
 861         if (cluster0->GetRow()>63) ipad=2;
 862       }
 863       //
 864       Float_t sigmaY =msigmay;
 865       Float_t sigmaZ =msigmay;
 866       if (iter==2){
 867         sigmaY = fClusterParam->GetErrorQParScaled(0,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 868                                                            TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
 869         sigmaZ = fClusterParam->GetErrorQParScaled(1,ipad,(250.0-TMath::Abs(cluster0->GetZ())),
 870                                                            TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
 871       }
 872       Double_t deltaX = cluster0->GetX()-refX;
 873       Double_t x[5];
 874       x[0] = (ipad==0) ? 0:1;
 875       x[1] = deltaX;
 876       x[2] = (ipad==0) ? 0:deltaX;
 877       x[3] = deltaX*deltaX;
 878       if (ipad<2){
 879         fitterY.AddPoint(x,y,sigmaY);
 880         fitterZ.AddPoint(x,z,sigmaZ);
 881       }
 882       if (ipad==0){
 883         nclDown++;
 884         fitterY0.AddPoint(&deltaX,y,sigmaY);
 885         fitterZ0.AddPoint(&deltaX,z,sigmaZ);
 886       }
 887       if (ipad==1){
 888         nclUp++;
 889         fitterY1.AddPoint(&deltaX,y,sigmaY);
 890         fitterZ1.AddPoint(&deltaX,z,sigmaZ);
 891       }
 892     }
 893     if (nclUp<kMinClustersSector) continue;
 894     if (nclDown<kMinClustersSector) continue;
 895     fitterY.Eval();
 896     fitterZ.Eval();
 897     fitterY0.Eval();
 898     fitterZ0.Eval();
 899     fitterY1.Eval();
 900     fitterZ1.Eval();
 901     param0[0] = fitterY0.GetParameter(0);
 902     param0[1] = fitterY0.GetParameter(1);
 903     param0[2] = fitterY0.GetParameter(2);
 904     param1[0] = fitterY1.GetParameter(0);
 905     param1[1] = fitterY1.GetParameter(1);
 906     param1[2] = fitterY1.GetParameter(2);
 907     //
 908     angley = fitterY.GetParameter(2);
 909     anglez = fitterZ.GetParameter(2);
 910     //
 911     TVectorD    parY(5);
 912     TMatrixD    matY(5,5);
 913     TVectorD    parZ(5);
 914     TMatrixD    matZ(5,5);
 915     Double_t    chi2Y= fitterY.GetChisquare()/(nclUp+nclDown);
 916     Double_t    chi2Z= fitterZ.GetChisquare()/(nclUp+nclDown);
 917     fitterY.GetParameters(parY);
 918     fitterY.GetCovarianceMatrix(matY);
 919     fitterZ.GetParameters(parZ);
 920     fitterZ.GetCovarianceMatrix(matZ);
 921     if (iter==0) {
 922       msigmay = msigmay*TMath::Sqrt(chi2Y);
 923       msigmaz = msigmaz*TMath::Sqrt(chi2Z);
 924     }
 925     Float_t sigmaY  = TMath::Sqrt(matY(1,1)*chi2Y);
 926     Float_t sigmaDY = TMath::Sqrt(matY(3,3)*chi2Y);
 927     Float_t sigmaDDY = TMath::Sqrt(matY(4,4)*chi2Y);
 928     Float_t sigmaZ  = TMath::Sqrt(matZ(1,1)*chi2Z);
 929     Float_t sigmaDZ = TMath::Sqrt(matZ(3,3)*chi2Z);
 930     Float_t sigmaDDZ = TMath::Sqrt(matZ(4,4)*chi2Z);
 931     //
 932     TVectorD    parY0(3);
 933     TMatrixD    matY0(3,3);
 934     TVectorD    parZ0(3);
 935     TMatrixD    matZ0(3,3);
 936     Double_t    chi2Y0= fitterY0.GetChisquare()/(nclDown);
 937     Double_t    chi2Z0= fitterZ0.GetChisquare()/(nclDown);
 938     fitterY0.GetParameters(parY0);
 939     fitterY0.GetCovarianceMatrix(matY0);
 940     fitterZ0.GetParameters(parZ0);
 941     fitterZ0.GetCovarianceMatrix(matZ0);
 942     Float_t sigmaY0  = TMath::Sqrt(matY0(0,0)*chi2Y0);
 943     Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1)*chi2Y0);
 944     Float_t sigmaDDY0 = TMath::Sqrt(matY0(2,2)*chi2Y0);
 945     Float_t sigmaZ0  = TMath::Sqrt(matZ0(0,0)*chi2Z0);
 946     Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1)*chi2Z0);
 947     Float_t sigmaDDZ0 = TMath::Sqrt(matZ0(2,2)*chi2Z0);
 948     //
 949     TVectorD    parY1(3);
 950     TMatrixD    matY1(3,3);
 951     TVectorD    parZ1(3);
 952     TMatrixD    matZ1(3,3);
 953     Double_t    chi2Y1= fitterY1.GetChisquare()/(nclUp);
 954     Double_t    chi2Z1= fitterZ1.GetChisquare()/(nclUp);
 955     fitterY1.GetParameters(parY1);
 956     fitterY1.GetCovarianceMatrix(matY1);
 957     fitterZ1.GetParameters(parZ1);
 958     fitterZ1.GetCovarianceMatrix(matZ1);
 959     Float_t sigmaY1  = TMath::Sqrt(matY1(0,0)*chi2Y1);
 960     Float_t sigmaDY1 = TMath::Sqrt(matY1(1,1)*chi2Y1);
 961     Float_t sigmaDDY1 = TMath::Sqrt(matY1(2,2)*chi2Y1);
 962     Float_t sigmaZ1  = TMath::Sqrt(matZ1(0,0)*chi2Z1);
 963     Float_t sigmaDZ1 = TMath::Sqrt(matZ1(1,1)*chi2Z1);
 964     Float_t sigmaDDZ1 = TMath::Sqrt(matZ1(2,2)*chi2Z1);
 965     const AliESDfriendTrack * ftrack = esdTrack->GetFriendTrack();
 966     AliTrackPointArray *points = (AliTrackPointArray*)ftrack->GetTrackPointArray();
 967
 968     if (iter>0) (*fDebugStream)<<"Align"<<
 969       "track.="<<track<<
 970       "Iter="<<iter<<
 971       "xref="<<refX<<
 972       "Points="<<points<<
 973       "Sector="<<rSector<<
 974       "nclUp="<<nclUp<<
 975       "nclDown="<<nclDown<<
 976       "angley="<<angley<<
 977       "anglez="<<anglez<<
 978       //
 979       "chi2Y="<<chi2Y<<
 980       "chi2Z="<<chi2Z<<
 981       "parY.="<<&parY<<
 982       "parZ.="<<&parZ<<
 983       "matY.="<<&matY<<
 984       "matZ.="<<&matZ<<
 985       "sigmaY="<<sigmaY<<
 986       "sigmaZ="<<sigmaZ<<
 987       "sigmaDY="<<sigmaDY<<
 988       "sigmaDZ="<<sigmaDZ<<
 989       "sigmaDDY="<<sigmaDDY<<
 990       "sigmaDDZ="<<sigmaDDZ<<
 991       //
 992       "chi2Y0="<<chi2Y0<<
 993       "chi2Z0="<<chi2Z0<<
 994       "parY0.="<<&parY0<<
 995       "parZ0.="<<&parZ0<<
 996       "matY0.="<<&matY0<<
 997       "matZ0.="<<&matZ0<<
 998       "sigmaY0="<<sigmaY0<<
 999       "sigmaZ0="<<sigmaZ0<<
1000       "sigmaDY0="<<sigmaDY0<<
1001       "sigmaDZ0="<<sigmaDZ0<<
1002       "sigmaDDY0="<<sigmaDDY0<<
1003       "sigmaDDZ0="<<sigmaDDZ0<<
1004       //
1005       "chi2Y1="<<chi2Y1<<
1006       "chi2Z1="<<chi2Z1<<
1007       "parY1.="<<&parY1<<
1008       "parZ1.="<<&parZ1<<
1009       "matY1.="<<&matY1<<
1010       "matZ1.="<<&matZ1<<
1011       "sigmaY1="<<sigmaY1<<
1012       "sigmaZ1="<<sigmaZ1<<
1013       "sigmaDY1="<<sigmaDY1<<
1014       "sigmaDZ1="<<sigmaDZ1<<
1015       "sigmaDDY1="<<sigmaDDY1<<
1016       "sigmaDDZ1="<<sigmaDDZ1<<
1017       "\n";
1018   }
1019 }