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