TPC/Calib/AliTPCcalibAlign.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
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //     Class to make a internal alignemnt of TPC chambers                    //
  20 //
  21 //     Requierements - Warnings:
  22 //     1. Before using this componenent the magnetic filed has to be set properly //
  23 //     2. The systematic effects  - unlinearities has to be understood
  24 //
  25 //     If systematic and unlinearities are not under control
  26 //     the alignment is just effective alignment. Not second order corrction
  27 //     are calculated.
  28 //
  29 //     The histograming of the edge effects and unlineratities integral part
  30 //     of the component (currently only in debug stream)
  31 //
  32 //     3 general type of linear transformation investigated (see bellow)
  33 //
  34 //     By default only 6 parameter alignment to be used - other just for QA purposes
  35
  36 //     Different linear tranformation investigated
  37 //     12 parameters - arbitrary linear transformation
  38 //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
  39 //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
  40 //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
  41 //
  42 //      9 parameters - scaling fixed to 1
  43 //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
  44 //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
  45 //                     a20  a21  a22 a23     p[4]   p[5]  1      p[8]
  46 //
  47 //      6 parameters - x-y rotation x-z, y-z tiliting
  48 //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
  49 //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
  50 //                     a20  a21  a22 a23     p[1]   p[2]  1     p[5]
  51 //
  52 //
  53 //      Debug stream supported
  54 //      0. Align    - The main output of the Alignment component
  55 //                  - Used for visualization of the misalignment between sectors
  56 //                  - Results of the missalignment fit and the mean and sigmas of histograms
  57 //                   stored there
  58 //      1. Tracklet - StreamLevel >1
  59 //                  - Dump all information about tracklet match from sector1 to sector 2
  60 //                  - Default histogram residulas created in parallel
  61 //                  - Check this streamer in case of suspicious content of these histograms
  62 //      2. Track    - StreamLevel>5
  63 //                  - For debugging of the edge effects
  64 //                  - All information  - extrapolation inside of one sectors
  65 //                  - Created in order to distinguish between unlinearities inside of o
  66 //                    sector and  missalignment
  67
  68 //
  69 //
  70 /*
  71   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  72   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  73   AliXRDPROOFtoolkit tool;
  74   TChain * chain = tool.MakeChain("align.txt","Track",0,10200);
  75   chain->Lookup();
  76   TCut cutA("abs(tp1.fP[1]-tp2.fP[1])<0.3&&abs(tp1.fP[0]-tp2.fP[0])<0.15&&abs(tp1.fP[3]-tp2.fP[3])<0.01&&abs(tp1.fP[2]-tp2.fP[2])<0.01");
  77   TCut cutS("s1%36==s2%36");
  78
  79   .x ~/UliStyle.C
  80   .x $ALICE_ROOT/macros/loadlibsREC.C
  81
  82   gSystem->Load("$ROOTSYS/lib/libXrdClient.so");
  83   gSystem->Load("libProof");
  84   gSystem->Load("libANALYSIS");
  85   gSystem->Load("libSTAT");
  86   gSystem->Load("libTPCcalib");
  87   //
  88   // compare reference
  89   TFile fcalib("CalibObjects.root");
  90   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  91
  92   AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
  93   //
  94   //
  95   align->EvalFitters();
  96   align->MakeTree("alignTree.root");
  97   TFile falignTree("alignTree.root");
  98   TTree * treeAlign = (TTree*)falignTree.Get("Align");
  99
 100
 101 */
 102
 103 ////
 104 ////
 105
 106 #include "TLinearFitter.h"
 107 #include "AliTPCcalibAlign.h"
 108 #include "AliTPCROC.h"
 109 #include "AliTPCPointCorrection.h"
 110 #include "AliTrackPointArray.h"
 111
 112 #include "AliExternalTrackParam.h"
 113 #include "AliESDEvent.h"
 114 #include "AliESDfriend.h"
 115 #include "AliESDtrack.h"
 116
 117 #include "AliTPCTracklet.h"
 118 #include "TH1D.h"
 119 #include "TH2F.h"
 120 #include "THnSparse.h"
 121 #include "THn.h"
 122 #include "TVectorD.h"
 123 #include "TTreeStream.h"
 124 #include "TFile.h"
 125 #include "TTree.h"
 126 #include "TF1.h"
 127 #include "TGraphErrors.h"
 128 #include "AliTPCclusterMI.h"
 129 #include "AliTPCseed.h"
 130 #include "AliTracker.h"
 131 #include "TClonesArray.h"
 132 #include "AliLog.h"
 133 #include "TFile.h"
 134 #include "TProfile.h"
 135 #include "TCanvas.h"
 136 #include "TDatabasePDG.h"
 137
 138 #include "TTreeStream.h"
 139 #include "Riostream.h"
 140 #include "TRandom.h"
 141 #include <sstream>
 142
 143 #include "AliSysInfo.h"
 144 using namespace std;
 145
 146 AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0;
 147 Double_t          AliTPCcalibAlign::fgkMergeEntriesCut=10000000.; //10**7 tracks
 148 ClassImp(AliTPCcalibAlign)
 149
 150
 151
 152
 153 AliTPCcalibAlign* AliTPCcalibAlign::Instance()
 154 {
 155   //
 156   // Singleton implementation
 157   // Returns an instance of this class, it is created if neccessary
 158   //
 159   if (fgInstance == 0){
 160     fgInstance = new AliTPCcalibAlign();
 161   }
 162   return fgInstance;
 163 }
 164
 165
 166
 167
 168 AliTPCcalibAlign::AliTPCcalibAlign()
 169   :  AliTPCcalibBase(),
 170      fDphiHistArray(72*72),
 171      fDthetaHistArray(72*72),
 172      fDyHistArray(72*72),
 173      fDzHistArray(72*72),
 174      //
 175      fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
 176      fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
 177      fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
 178      fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
 179      fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
 180      fDzZHistArray(72*72),        // array of residual histograms  z     -kZz
 181      fFitterArray12(72*72),
 182      fFitterArray9(72*72),
 183      fFitterArray6(72*72),
 184      fMatrixArray12(72*72),
 185      fMatrixArray9(72*72),
 186      fMatrixArray6(72*72),
 187      fCombinedMatrixArray6(72),
 188      fNoField(kFALSE),
 189      fXIO(0),
 190      fXmiddle(0),
 191      fXquadrant(0),
 192      fArraySectorIntParam(36), // array of sector alignment parameters
 193      fArraySectorIntCovar(36), // array of sector alignment covariances
 194      //
 195      // Kalman filter for global alignment
 196      //
 197      fSectorParamA(0),     // Kalman parameter   for A side
 198      fSectorCovarA(0),     // Kalman covariance  for A side
 199      fSectorParamC(0),     // Kalman parameter   for A side
 200      fSectorCovarC(0),     // Kalman covariance  for A side
 201      fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
 202 {
 203   //
 204   // Constructor
 205   //
 206   for (Int_t i=0;i<72*72;++i) {
 207     fPoints[i]=0;
 208   }
 209   AliTPCROC * roc = AliTPCROC::Instance();
 210   fXquadrant = roc->GetPadRowRadii(36,53);
 211   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
 212   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 213   fClusterDelta[0]=0;   // cluster residuals -  Y
 214   fClusterDelta[1]=0;   // cluster residuals -  Z
 215
 216
 217   fTrackletDelta[0]=0;   // tracklet residuals
 218   fTrackletDelta[1]=0;   // tracklet residuals
 219   fTrackletDelta[2]=0;   // tracklet residuals
 220   fTrackletDelta[3]=0;   // tracklet residuals
 221 }
 222
 223 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
 224   :AliTPCcalibBase(),
 225    fDphiHistArray(72*72),
 226    fDthetaHistArray(72*72),
 227    fDyHistArray(72*72),
 228    fDzHistArray(72*72),
 229    fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
 230    fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
 231    fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
 232    fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
 233    fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
 234    fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     //
 235    fFitterArray12(72*72),
 236    fFitterArray9(72*72),
 237    fFitterArray6(72*72),
 238    fMatrixArray12(72*72),
 239    fMatrixArray9(72*72),
 240    fMatrixArray6(72*72),
 241    fCombinedMatrixArray6(72),
 242    fNoField(kFALSE),
 243    fXIO(0),
 244    fXmiddle(0),
 245    fXquadrant(0),
 246    fArraySectorIntParam(36), // array of sector alignment parameters
 247    fArraySectorIntCovar(36), // array of sector alignment covariances
 248    //
 249    // Kalman filter for global alignment
 250    //
 251    fSectorParamA(0),     // Kalman parameter   for A side
 252    fSectorCovarA(0),     // Kalman covariance  for A side
 253    fSectorParamC(0),     // Kalman parameter   for A side
 254    fSectorCovarC(0),     // Kalman covariance  for A side
 255    fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
 256
 257 {
 258   //
 259   // Constructor
 260   //
 261   SetName(name);
 262   SetTitle(title);
 263   for (Int_t i=0;i<72*72;++i) {
 264     fPoints[i]=0;
 265   }
 266   AliTPCROC * roc = AliTPCROC::Instance();
 267   fXquadrant = roc->GetPadRowRadii(36,53);
 268   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
 269   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 270   fClusterDelta[0]=0;   // cluster residuals
 271   fClusterDelta[1]=0;   // cluster residuals
 272
 273   fTrackletDelta[0]=0;   // tracklet residuals
 274   fTrackletDelta[1]=0;   // tracklet residuals
 275   fTrackletDelta[2]=0;   // tracklet residuals
 276   fTrackletDelta[3]=0;   // tracklet residuals
 277 }
 278
 279
 280 AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
 281   :AliTPCcalibBase(align),
 282    fDphiHistArray(align.fDphiHistArray),
 283    fDthetaHistArray(align.fDthetaHistArray),
 284    fDyHistArray(align.fDyHistArray),
 285    fDzHistArray(align.fDzHistArray),
 286    fDyPhiHistArray(align.fDyPhiHistArray),      // array of residual histograms  y     -kYPhi
 287    fDzThetaHistArray(align.fDzThetaHistArray),    // array of residual histograms  z-z   -kZTheta
 288    fDphiZHistArray(align.fDphiZHistArray),      // array of residual histograms  phi   -kPhiz
 289    fDthetaZHistArray(align.fDthetaZHistArray),    // array of residual histograms  theta -kThetaz
 290    fDyZHistArray(align.fDyZHistArray),        // array of residual histograms  y     -kYz
 291    fDzZHistArray(align.fDzZHistArray),        // array of residual histograms  z     -kZz
 292    //
 293    fFitterArray12(align.fFitterArray12),
 294    fFitterArray9(align.fFitterArray9),
 295    fFitterArray6(align.fFitterArray6),
 296
 297    fMatrixArray12(align.fMatrixArray12),
 298    fMatrixArray9(align.fMatrixArray9),
 299    fMatrixArray6(align.fMatrixArray6),
 300    fCombinedMatrixArray6(align.fCombinedMatrixArray6),
 301    fNoField(align.fNoField),
 302    fXIO(align.fXIO),
 303    fXmiddle(align.fXmiddle),
 304    fXquadrant(align.fXquadrant),
 305    fArraySectorIntParam(align.fArraySectorIntParam), // array of sector alignment parameters
 306    fArraySectorIntCovar(align.fArraySectorIntCovar), // array of sector alignment covariances
 307    fSectorParamA(0),     // Kalman parameter   for A side
 308    fSectorCovarA(0),     // Kalman covariance  for A side
 309    fSectorParamC(0),     // Kalman parameter   for A side
 310    fSectorCovarC(0),      // Kalman covariance  for A side
 311   fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
 312
 313 {
 314   //
 315   // copy constructor - copy also the content
 316   //
 317   TH1 * his = 0;
 318   TObjArray * arr0=0;
 319   const TObjArray *arr1=0;
 320   for (Int_t index =0; index<72*72; index++){
 321     for (Int_t iarray=0;iarray<10; iarray++){
 322       if (iarray==kY){
 323         arr0 = &fDyHistArray;
 324         arr1 = &align.fDyHistArray;
 325       }
 326       if (iarray==kZ){
 327         arr0 = &fDzHistArray;
 328         arr1 = &align.fDzHistArray;
 329       }
 330       if (iarray==kPhi){
 331         arr0 = &fDphiHistArray;
 332         arr1 = &align.fDphiHistArray;
 333       }
 334       if (iarray==kTheta){
 335         arr0 = &fDthetaHistArray;
 336         arr1 = &align.fDthetaHistArray;
 337       }
 338       if (iarray==kYz){
 339         arr0 = &fDyZHistArray;
 340         arr1 = &align.fDyZHistArray;
 341       }
 342       if (iarray==kZz){
 343         arr0 = &fDzZHistArray;
 344         arr1 = &align.fDzZHistArray;
 345       }
 346       if (iarray==kPhiZ){
 347         arr0 = &fDphiZHistArray;
 348         arr1 = &align.fDphiZHistArray;
 349       }
 350       if (iarray==kThetaZ){
 351         arr0 = &fDthetaZHistArray;
 352         arr1 = &align.fDthetaZHistArray;
 353       }
 354
 355       if (iarray==kYPhi){
 356         arr0 = &fDyPhiHistArray;
 357         arr1 = &align.fDyPhiHistArray;
 358       }
 359       if (iarray==kZTheta){
 360         arr0 = &fDzThetaHistArray;
 361         arr1 = &align.fDzThetaHistArray;
 362       }
 363
 364       if (arr1->At(index)) {
 365         his = (TH1*)arr1->At(index)->Clone();
 366         his->SetDirectory(0);
 367         arr0->AddAt(his,index);
 368       }
 369     }
 370   }
 371   //
 372   //
 373   //
 374   if (align.fSectorParamA){
 375     fSectorParamA = (TMatrixD*)align.fSectorParamA->Clone();
 376     fSectorParamA = (TMatrixD*)align.fSectorCovarA->Clone();
 377     fSectorParamC = (TMatrixD*)align.fSectorParamA->Clone();
 378     fSectorParamC = (TMatrixD*)align.fSectorCovarA->Clone();
 379   }
 380   fClusterDelta[0]=0;   // cluster residuals
 381   fClusterDelta[1]=0;   // cluster residuals
 382
 383   fTrackletDelta[0]=0;   // tracklet residuals
 384   fTrackletDelta[1]=0;   // tracklet residuals
 385   fTrackletDelta[2]=0;   // tracklet residuals
 386   fTrackletDelta[3]=0;   // tracklet residuals
 387 }
 388
 389
 390 AliTPCcalibAlign::~AliTPCcalibAlign() {
 391   //
 392   // destructor
 393   //
 394   fDphiHistArray.SetOwner(kTRUE);    // array of residual histograms  phi      -kPhi
 395   fDthetaHistArray.SetOwner(kTRUE);  // array of residual histograms  theta    -kTheta
 396   fDyHistArray.SetOwner(kTRUE);      // array of residual histograms  y        -kY
 397   fDzHistArray.SetOwner(kTRUE);      // array of residual histograms  z        -kZ
 398   //
 399   fDyPhiHistArray.SetOwner(kTRUE);      // array of residual histograms  y     -kYPhi
 400   fDzThetaHistArray.SetOwner(kTRUE);    // array of residual histograms  z-z   -kZTheta
 401   //
 402   fDphiZHistArray.SetOwner(kTRUE);      // array of residual histograms  phi   -kPhiz
 403   fDthetaZHistArray.SetOwner(kTRUE);    // array of residual histograms  theta -kThetaz
 404   fDyZHistArray.SetOwner(kTRUE);        // array of residual histograms  y     -kYz
 405   fDzZHistArray.SetOwner(kTRUE);        // array of residual histograms  z     -kZz
 406
 407   fDphiHistArray.Delete();    // array of residual histograms  phi      -kPhi
 408   fDthetaHistArray.Delete();  // array of residual histograms  theta    -kTheta
 409   fDyHistArray.Delete();      // array of residual histograms  y        -kY
 410   fDzHistArray.Delete();      // array of residual histograms  z        -kZ
 411   //
 412   fDyPhiHistArray.Delete();      // array of residual histograms  y     -kYPhi
 413   fDzThetaHistArray.Delete();    // array of residual histograms  z-z   -kZTheta
 414   //
 415   fDphiZHistArray.Delete();      // array of residual histograms  phi   -kPhiz
 416   fDthetaZHistArray.Delete();    // array of residual histograms  theta -kThetaz
 417   fDyZHistArray.Delete();        // array of residual histograms  y     -kYz
 418   fDzZHistArray.Delete();        // array of residual histograms  z     -kZz
 419
 420   fFitterArray12.SetOwner(kTRUE);    // array of fitters
 421   fFitterArray9.SetOwner(kTRUE);     // array of fitters
 422   fFitterArray6.SetOwner(kTRUE);     // array of fitters
 423   //
 424   fMatrixArray12.SetOwner(kTRUE);    // array of transnformtation matrix
 425   fMatrixArray9.SetOwner(kTRUE);     // array of transnformtation matrix
 426   fMatrixArray6.SetOwner(kTRUE);     // array of transnformtation matrix
 427   //
 428   fFitterArray12.Delete();    // array of fitters
 429   fFitterArray9.Delete();     // array of fitters
 430   fFitterArray6.Delete();     // array of fitters
 431   //
 432   fMatrixArray12.Delete();    // array of transnformtation matrix
 433   fMatrixArray9.Delete();     // array of transnformtation matrix
 434   fMatrixArray6.Delete();     // array of transnformtation matrix
 435
 436
 437   fArraySectorIntParam.SetOwner(kTRUE); // array of sector alignment parameters
 438   fArraySectorIntCovar.SetOwner(kTRUE); // array of sector alignment covariances
 439   fArraySectorIntParam.Delete(); // array of sector alignment parameters
 440   fArraySectorIntCovar.Delete(); // array of sector alignment covariances
 441   for (Int_t i=0; i<2; i++){
 442     delete fClusterDelta[i];   // cluster residuals
 443   }
 444
 445   for (Int_t i=0; i<4; i++){
 446     delete fTrackletDelta[i];   // tracklet residuals
 447   }
 448
 449
 450 }
 451
 452 void AliTPCcalibAlign::Process(AliESDEvent *event) {
 453   //
 454   // Process pairs of cosmic tracks
 455   //
 456   const Double_t kptDownscale=50; // downscale factor for the low pt particels
 457   if (!fClusterDelta[0])  MakeResidualHistos();
 458   if (!fTrackletDelta[0])  MakeResidualHistosTracklet();
 459   //
 460   fCurrentEvent=event;
 461   ExportTrackPoints(event);  // export track points for external calibration
 462   const Int_t kMaxTracks =6;
 463   const Int_t kminCl = 40;
 464   AliESDfriend *eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 465   if (!eESDfriend) return;
 466   if (eESDfriend->TestSkipBit()) return;
 467   Int_t ntracks=event->GetNumberOfTracks();
 468   Float_t dca0[2];
 469   Float_t dca1[2];
 470   //
 471   //
 472   // process seeds
 473   //
 474   for (Int_t i0=0;i0<ntracks;++i0) {
 475     AliESDtrack *track0 = event->GetTrack(i0);
 476     AliESDfriendTrack *friendTrack = 0;
 477     TObject *calibObject=0;
 478     AliTPCseed *seed0 = 0;
 479     //
 480     friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
 481     if (!friendTrack) continue;
 482     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 483       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
 484     }
 485     if (!seed0) continue;
 486     fCurrentTrack=track0;
 487     fCurrentFriendTrack=friendTrack;
 488     fCurrentSeed=seed0;
 489     fCurrentEvent=event;
 490     Double_t scalept= TMath::Min(1./TMath::Abs(track0->GetParameter()[4]),2.);
 491     Bool_t   isSelected = (TMath::Exp(2*scalept)>kptDownscale*gRandom->Rndm());
 492     if (isSelected) ProcessSeed(seed0);
 493   }
 494   //
 495   // process cosmic pairs
 496   //
 497   if (ntracks>kMaxTracks) return;
 498   //
 499   //select pairs - for alignment
 500   for (Int_t i0=0;i0<ntracks;++i0) {
 501     AliESDtrack *track0 = event->GetTrack(i0);
 502     //    if (track0->GetTPCNcls()<kminCl) continue;
 503     track0->GetImpactParameters(dca0[0],dca0[1]);
 504     //    if (TMath::Abs(dca0[0])>30) continue;
 505     //
 506     for (Int_t i1=0;i1<ntracks;++i1) {
 507       if (i0==i1) continue;
 508       AliESDtrack *track1 = event->GetTrack(i1);
 509       //      if (track1->GetTPCNcls()<kminCl) continue;
 510       track1->GetImpactParameters(dca1[0],dca1[1]);
 511       // fast cuts on dca and theta
 512       //      if (TMath::Abs(dca1[0]+dca0[0])>15) continue;
 513       //      if (TMath::Abs(dca1[1]-dca0[1])>15) continue;
 514       if (TMath::Abs(track0->GetParameter()[3]+track1->GetParameter()[3])>0.1) continue;
 515       //
 516       AliESDfriendTrack *friendTrack = 0;
 517       TObject *calibObject=0;
 518       AliTPCseed *seed0 = 0,*seed1=0;
 519       //
 520       friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
 521       if (!friendTrack) continue;
 522       for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 523         if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
 524       }
 525       friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i1);;
 526       if (!friendTrack) continue;
 527       for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 528         if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
 529       }
 530       if (!seed0) continue;
 531       //
 532       //
 533       if (!seed1) continue;
 534       Int_t nclsectors0[72], nclsectors1[72];
 535       for (Int_t isec=0;isec<72;isec++){
 536         nclsectors0[isec]=0;
 537         nclsectors1[isec]=0;
 538       }
 539       for (Int_t i=0;i<160;i++){
 540         AliTPCclusterMI *c0=seed0->GetClusterPointer(i);
 541         AliTPCclusterMI *c1=seed1->GetClusterPointer(i);
 542         if (c0)  nclsectors0[c0->GetDetector()]+=1;
 543         if (c1)  nclsectors1[c1->GetDetector()]+=1;
 544       }
 545
 546       for (Int_t isec0=0; isec0<72;isec0++){
 547         if (nclsectors0[isec0]<kminCl) continue;
 548         for (Int_t isec1=0; isec1<72;isec1++){
 549           if (nclsectors1[isec1]<kminCl) continue;
 550           Int_t s0 = isec0;
 551           Int_t s1 = isec1;
 552           Double_t parLine0[10];
 553           Double_t parLine1[10];
 554           TMatrixD par0(4,1),cov0(4,4),par1(4,1),cov1(4,4);
 555           Bool_t useInnerOuter = kFALSE;
 556           if (s1%36!=s0%36) useInnerOuter = fUseInnerOuter;  // for left - right alignment both sectors refit can be used if specified
 557           Int_t nl0 = RefitLinear(seed0,s0, parLine0, s0,par0,cov0,fXIO,useInnerOuter);
 558           Int_t nl1 = RefitLinear(seed1,s1, parLine1, s0,par1,cov1,fXIO,useInnerOuter);
 559           parLine0[0]=0;  // reference frame in IO boundary
 560           parLine1[0]=0;
 561           //      if (nl0<kminCl || nl1<kminCl) continue;
 562           //
 563           //
 564           Bool_t isOK=kTRUE;
 565           if (TMath::Min(nl0,nl1)<kminCl) isOK=kFALSE;
 566           // apply selection criteria
 567           //
 568           Float_t dp0,dp1,dp3;
 569           Float_t pp0,pp1,pp3;
 570           dp0=par0(0,0)-par1(0,0);
 571           dp1=par0(1,0)-par1(1,0);
 572           dp3=par0(3,0)-par1(3,0);
 573           pp0=dp0/TMath::Sqrt(cov0(0,0)+cov1(0,0)+0.1*0.1);
 574           pp1=dp1/TMath::Sqrt(cov0(1,1)+cov1(1,1)+0.0015*0.0015);
 575           pp3=dp3/TMath::Sqrt(cov0(3,3)+cov1(3,3)+0.0015*0.0015);
 576           //
 577           if (TMath::Abs(dp0)>1.0)  isOK=kFALSE;
 578           if (TMath::Abs(dp1)>0.02) isOK=kFALSE;
 579           if (TMath::Abs(dp3)>0.02) isOK=kFALSE;
 580           if (TMath::Abs(pp0)>6)  isOK=kFALSE;
 581           if (TMath::Abs(pp1)>6) isOK=kFALSE;
 582           if (TMath::Abs(pp3)>6) isOK=kFALSE;
 583           //
 584           if (isOK){
 585             FillHisto(parLine0,parLine1,s0,s1);
 586             ProcessAlign(parLine0,parLine1,s0,s1);
 587             UpdateKalman(s0,s1,par0, cov0, par1, cov1);
 588           }
 589           if (fStreamLevel>0){
 590             TTreeSRedirector *cstream = GetDebugStreamer();
 591             if (cstream){
 592               (*cstream)<<"cosmic"<<
 593                 "isOK="<<isOK<<
 594                 "s0="<<s0<<
 595                 "s1="<<s1<<
 596                 "nl0="<<nl0<<
 597                 "nl1="<<nl1<<
 598                 "p0.="<<&par0<<
 599                 "p1.="<<&par1<<
 600                 "c0.="<<&cov0<<
 601                 "c1.="<<&cov1<<
 602                 "\n";
 603             }
 604           }
 605         }
 606       }
 607     }
 608   }
 609 }
 610
 611 void  AliTPCcalibAlign::ExportTrackPoints(AliESDEvent *event){
 612   //
 613   // Export track points for alignment - calibration
 614   // export space points for pairs of tracks if possible
 615   //
 616   AliESDfriend *eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 617   if (!eESDfriend) return;
 618   Int_t ntracks=event->GetNumberOfTracks();
 619   Int_t kMaxTracks=4;   // maximal number of tracks for cosmic pairs
 620   Int_t kMinVertexTracks=5;   // maximal number of tracks for vertex mesurement
 621
 622   //cuts
 623   const Int_t kminCl     = 60;
 624   const Int_t kminClSum  = 120;
 625   //const Double_t kDistY  = 5;
 626   // const Double_t kDistZ  = 40;
 627   const Double_t kDistTh = 0.05;
 628   const Double_t kDistThS = 0.002;
 629   const Double_t kDist1Pt = 0.1;
 630   const Double_t kMaxD0 =3;  // max distance to the primary vertex
 631   const Double_t kMaxD1 =5;  // max distance to the primary vertex
 632   const AliESDVertex *tpcVertex = 0;
 633   // get the primary vertex TPC
 634   if (ntracks>kMinVertexTracks) {
 635     tpcVertex = event->GetPrimaryVertexSPD();
 636     if (tpcVertex->GetNContributors()<kMinVertexTracks) tpcVertex=0;
 637   }
 638   //
 639   Float_t dca0[2];
 640   //  Float_t dca1[2];
 641   Int_t index0=0,index1=0;
 642   //
 643   for (Int_t i0=0;i0<ntracks;++i0) {
 644     AliESDtrack *track0 = event->GetTrack(i0);
 645     if (!track0) continue;
 646     if ((track0->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
 647     if (track0->GetOuterParam()==0) continue;
 648     if (track0->GetInnerParam()==0) continue;
 649     if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt()-track0->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
 650     if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
 651     if (TMath::Abs(track0->GetInnerParam()->GetTgl()-track0->GetOuterParam()->GetTgl())>kDistThS) continue;
 652     AliESDtrack *track1P = 0;
 653     if (track0->GetTPCNcls()<kminCl) continue;
 654     track0->GetImpactParameters(dca0[0],dca0[1]);
 655     index0=i0;
 656     index1=-1;
 657     //
 658     if (ntracks<kMaxTracks) for (Int_t i1=i0+1;i1<ntracks;++i1) {
 659       if (i0==i1) continue;
 660       AliESDtrack *track1 = event->GetTrack(i1);
 661       if (!track1) continue;
 662       if ((track1->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
 663       if (track1->GetOuterParam()==0) continue;
 664       if (track1->GetInnerParam()==0) continue;
 665       if (track1->GetTPCNcls()<kminCl) continue;
 666       if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt()-track1->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
 667       if (TMath::Abs(track1->GetInnerParam()->GetTgl()-track1->GetOuterParam()->GetTgl())>kDistThS) continue;
 668       if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
 669       //track1->GetImpactParameters(dca1[0],dca1[1]);
 670       //if (TMath::Abs(dca1[0]-dca0[0])>kDistY) continue;
 671       //if (TMath::Abs(dca1[1]-dca0[1])>kDistZ) continue;
 672       if (TMath::Abs(track0->GetTgl()+track1->GetTgl())>kDistTh) continue;
 673       if (TMath::Abs(track0->GetSigned1Pt()+track1->GetSigned1Pt())>kDist1Pt) continue;
 674       track1P = track1;
 675       index1=i1;
 676     }
 677     AliESDfriendTrack *friendTrack = 0;
 678     TObject *calibObject=0;
 679     AliTPCseed *seed0 = 0,*seed1=0;
 680     //
 681     friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index0);;
 682     if (!friendTrack) continue;
 683     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 684       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
 685     }
 686     if (index1>0){
 687       friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index1);;
 688       if (!friendTrack) continue;
 689       for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 690         if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
 691       }
 692     }
 693     //
 694     Int_t npoints=0, ncont=0;
 695     if (seed0) {npoints+=seed0->GetNumberOfClusters(); ncont++;}
 696     if (seed1) {npoints+=seed1->GetNumberOfClusters(); ncont++;}
 697     if (npoints<kminClSum) continue;
 698     Int_t cpoint=0;
 699     AliTrackPointArray array(npoints);
 700     if (tpcVertex){
 701       Double_t dxyz[3]={0,0,0};
 702       Double_t dc[6]={0,0,0};
 703       tpcVertex->GetXYZ(dxyz);
 704       tpcVertex->GetCovarianceMatrix(dc);
 705       Float_t xyz[3]={static_cast<Float_t>(dxyz[0]),static_cast<Float_t>(dxyz[1]),static_cast<Float_t>(dxyz[2])};
 706       Float_t cov[6]={static_cast<Float_t>(dc[0]+1),static_cast<Float_t>(dc[1]),static_cast<Float_t>(dc[2]+1),static_cast<Float_t>(dc[3]),static_cast<Float_t>(dc[4]),static_cast<Float_t>(dc[5]+100.)};
 707       AliTrackPoint point(xyz,cov,73);  // add point to not existing volume
 708       Float_t dtpc[2],dcov[3];
 709       track0->GetImpactParametersTPC(dtpc,dcov);
 710       if (TMath::Abs(dtpc[0])>kMaxD0) continue;
 711       if (TMath::Abs(dtpc[1])>kMaxD1) continue;
 712     }
 713
 714     if (seed0) for (Int_t icl = 0; icl<160; icl++){
 715       AliTPCclusterMI *cluster=seed0->GetClusterPointer(icl);
 716       if (!cluster) continue;
 717       Float_t xyz[3];
 718       Float_t cov[6];
 719       cluster->GetGlobalXYZ(xyz);
 720       cluster->GetGlobalCov(cov);
 721       AliTrackPoint point(xyz,cov,cluster->GetDetector());
 722       array.AddPoint(npoints, &point);
 723       if (cpoint>=npoints) continue;  //shoul not happen
 724       array.AddPoint(cpoint, &point);
 725       cpoint++;
 726     }
 727     if (seed1) for (Int_t icl = 0; icl<160; icl++){
 728       AliTPCclusterMI *cluster=seed1->GetClusterPointer(icl);
 729       if (!cluster) continue;
 730       Float_t xyz[3];
 731       Float_t cov[6];
 732       cluster->GetGlobalXYZ(xyz);
 733       cluster->GetGlobalCov(cov);
 734       AliTrackPoint point(xyz,cov,cluster->GetDetector());
 735       array.AddPoint(npoints, &point);
 736       if (cpoint>=npoints) continue;  //shoul not happen
 737       array.AddPoint(cpoint, &point);
 738       cpoint++;
 739     }
 740     //
 741     //
 742     //
 743     TTreeSRedirector *cstream = GetDebugStreamer();
 744     if (cstream){
 745       Bool_t isVertex=(tpcVertex)? kTRUE:kFALSE;
 746       Double_t tof0=track0->GetTOFsignal();
 747       Double_t tof1=(track1P) ?  track1P->GetTOFsignal(): 0;
 748       static AliExternalTrackParam dummy;
 749       AliExternalTrackParam *p0In  = &dummy;
 750       AliExternalTrackParam *p1In  = &dummy;
 751       AliExternalTrackParam *p0Out = &dummy;
 752       AliExternalTrackParam *p1Out = &dummy;
 753       AliESDVertex vdummy;
 754       AliESDVertex *pvertex= (tpcVertex)? (AliESDVertex *)tpcVertex: &vdummy;
 755       if (track0) {
 756         p0In= new AliExternalTrackParam(*track0);
 757         p0Out=new AliExternalTrackParam(*(track0->GetOuterParam()));
 758       }
 759       if (track1P) {
 760         p1In= new AliExternalTrackParam(*track1P);
 761         p1Out=new AliExternalTrackParam(*(track1P->GetOuterParam()));
 762       }
 763
 764       (*cstream)<<"trackPoints"<<
 765         "run="<<fRun<<              //  run number
 766         "event="<<fEvent<<          //  event number
 767         "time="<<fTime<<            //  time stamp of event
 768         "trigger="<<fTrigger<<      //  trigger
 769         "triggerClass="<<&fTriggerClass<<      //  trigger
 770         "mag="<<fMagF<<             //  magnetic field
 771         "pvertex.="<<pvertex<<      // vertex
 772         //
 773         "isVertex="<<isVertex<<     // flag is with prim vertex
 774         "tof0="<<tof0<<             // tof signal 0
 775         "tof1="<<tof1<<             // tof signal 1
 776         "seed0.="<<seed0<<          // track info
 777         "ntracks="<<ntracks<<       // number of tracks
 778         "ncont="<<ncont<<           // number of contributors
 779         "p0In.="<<p0In<<              // track parameters0
 780         "p1In.="<<p1In<<              // track parameters1
 781         "p0Out.="<<p0Out<<          // track parameters0
 782         "p1Out.="<<p1Out<<          // track parameters0
 783         "p.="<<&array<<
 784         "\n";
 785     }
 786   }
 787 }
 788
 789
 790
 791
 792 void AliTPCcalibAlign::ProcessSeed(AliTPCseed *seed) {
 793   //
 794   //
 795   //
 796   // make a kalman tracklets out of seed
 797   //
 798   UpdateClusterDeltaField(seed);
 799   TObjArray tracklets=
 800     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
 801                                     kFALSE,20,4);
 802   tracklets.SetOwner();
 803   Int_t ntracklets = tracklets.GetEntries();
 804   if (ntracklets<2) return;
 805   //
 806   //
 807   for (Int_t i1=0;i1<ntracklets;i1++)
 808     for (Int_t i2=0;i2<ntracklets;i2++){
 809       if (i1==i2) continue;
 810       AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[i1]);
 811       AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[i2]);
 812       AliExternalTrackParam *common1=0,*common2=0;
 813       if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2)){
 814         ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
 815         UpdateAlignSector(seed,t1->GetSector());
 816       }
 817       delete common1;
 818       delete common2;
 819     }
 820 }
 821
 822 void AliTPCcalibAlign::Analyze(){
 823   //
 824   // Analyze function
 825   //
 826   EvalFitters();
 827 }
 828
 829
 830 void AliTPCcalibAlign::Terminate(){
 831   //
 832   // Terminate function
 833   // call base terminate + Eval of fitters
 834   //
 835   Info("AliTPCcalibAlign","Terminate");
 836   EvalFitters();
 837   AliTPCcalibBase::Terminate();
 838 }
 839
 840
 841 void AliTPCcalibAlign::UpdatePointCorrection(AliTPCPointCorrection * correction){
 842   //
 843   // Update point correction with alignment coefficients
 844   //
 845   for (Int_t isec=0;isec<36;isec++){
 846     TMatrixD * matCorr = (TMatrixD*)(correction->fArraySectorIntParam.At(isec));
 847     TMatrixD * matAlign = (TMatrixD*)(fArraySectorIntParam.At(isec));
 848     TMatrixD * matAlignCovar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
 849     if (!matAlign) continue;
 850     if (!matCorr) {
 851       correction->fArraySectorIntParam.AddAt(matAlign->Clone(),isec);
 852       correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec);
 853       continue;
 854     }
 855     (*matCorr)+=(*matAlign);
 856     correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec);
 857   }
 858   //
 859
 860 }
 861
 862
 863 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
 864                                         const AliExternalTrackParam &tp2,
 865                                         const AliTPCseed * seed,
 866                                         Int_t s1,Int_t s2) {
 867   //
 868   // Process function to fill fitters
 869   //
 870   if (!seed) return;
 871   Double_t t1[10],t2[10];
 872   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
 873   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
 874   x1   =tp1.GetX();
 875   y1   =tp1.GetY();
 876   z1   =tp1.GetZ();
 877   Double_t snp1=tp1.GetSnp();
 878   dydx1=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
 879   Double_t tgl1=tp1.GetTgl();
 880   // dz/dx = 1/(cos(theta)*cos(phi))
 881   dzdx1=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
 882   x2   =tp2.GetX();
 883   y2   =tp2.GetY();
 884   z2   =tp2.GetZ();
 885   Double_t snp2=tp2.GetSnp();
 886   dydx2=snp2/TMath::Sqrt((1.-snp2)*(1.+snp2));
 887   Double_t tgl2=tp2.GetTgl();
 888   dzdx2=tgl2/TMath::Sqrt((1.-snp2)*(1.+snp2));
 889
 890   //
 891   // Kalman parameters
 892   //
 893   t1[0]-=fXIO;
 894   t2[0]-=fXIO;
 895   // errors
 896   t1[5]=0; t2[5]=0;
 897   t1[6]=TMath::Sqrt(tp1.GetSigmaY2());
 898   t1[7]=TMath::Sqrt(tp1.GetSigmaSnp2());
 899   t1[8]=TMath::Sqrt(tp1.GetSigmaZ2());
 900   t1[9]=TMath::Sqrt(tp1.GetSigmaTgl2());
 901
 902   t2[6]=TMath::Sqrt(tp2.GetSigmaY2());
 903   t2[7]=TMath::Sqrt(tp2.GetSigmaSnp2());
 904   t2[8]=TMath::Sqrt(tp2.GetSigmaZ2());
 905   t2[9]=TMath::Sqrt(tp2.GetSigmaTgl2());
 906   //
 907   // linear parameters
 908   //
 909   Double_t parLine1[10]={0};
 910   Double_t parLine2[10]={0};
 911   TMatrixD par1(4,1),cov1(4,4),par2(4,1),cov2(4,4);
 912   Bool_t useInnerOuter = kFALSE;
 913   if (s1%36!=s2%36) useInnerOuter = fUseInnerOuter;  // for left - right alignment bot sectors refit can be used if specified
 914   Int_t nl1 = RefitLinear(seed,s1, parLine1, s1,par1,cov1,tp1.GetX(), useInnerOuter);
 915   Int_t nl2 = RefitLinear(seed,s2, parLine2, s1,par2,cov2,tp1.GetX(), useInnerOuter);
 916   parLine1[0]=tp1.GetX()-fXIO;   // parameters in  IROC-OROC boundary
 917   parLine2[0]=tp1.GetX()-fXIO;   // parameters in  IROC-OROC boundary
 918   //
 919   //
 920   //
 921   Int_t accept       =   AcceptTracklet(tp1,tp2);
 922   Int_t acceptLinear =   AcceptTracklet(parLine1,parLine2);
 923
 924
 925   if (fStreamLevel>1){
 926     TTreeSRedirector *cstream = GetDebugStreamer();
 927     if (cstream){
 928       static TVectorD vec1(5);
 929       static TVectorD vec2(5);
 930       static TVectorD vecL1(9);
 931       static TVectorD vecL2(9);
 932       vec1.SetElements(t1);
 933       vec2.SetElements(t2);
 934       vecL1.SetElements(parLine1);
 935       vecL2.SetElements(parLine2);
 936       AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
 937       AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
 938       (*cstream)<<"Tracklet"<<
 939         "accept="<<accept<<
 940         "acceptLinear="<<acceptLinear<<  // accept linear tracklets
 941         "run="<<fRun<<              //  run number
 942         "event="<<fEvent<<          //  event number
 943         "time="<<fTime<<            //  time stamp of event
 944         "trigger="<<fTrigger<<      //  trigger
 945         "triggerClass="<<&fTriggerClass<<      //  trigger
 946         "mag="<<fMagF<<             //  magnetic field
 947         "isOK="<<accept<<           //  flag - used for alignment
 948         "tp1.="<<p1<<
 949         "tp2.="<<p2<<
 950         "v1.="<<&vec1<<
 951         "v2.="<<&vec2<<
 952         "s1="<<s1<<
 953         "s2="<<s2<<
 954         "nl1="<<nl1<<       // linear fit - n points
 955         "nl2="<<nl2<<       // linear fit - n points
 956         "vl1.="<<&vecL1<<   // linear fits
 957         "vl2.="<<&vecL2<<   // linear fits
 958         "\n";
 959     }
 960   }
 961   if (TMath::Abs(fMagF)<0.005){
 962     //
 963     // use Linear fit
 964     //
 965     if (nl1>10 && nl2>10 &&(acceptLinear==0)){
 966       ProcessDiff(tp1,tp2, seed,s1,s2);
 967       if (TMath::Abs(parLine1[2])<0.8 &&TMath::Abs(parLine1[2])<0.8 ){ //angular cut
 968         FillHisto(parLine1,parLine2,s1,s2);
 969         ProcessAlign(parLine1,parLine2,s1,s2);
 970         FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
 971         FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
 972         //UpdateKalman(s1,s2,par1, cov1, par2, cov2); - OBSOLETE to be removed - 50 % of time here
 973       }
 974     }
 975   }
 976   if (accept>0) return;
 977   //
 978   // fill resolution histograms - previous cut included
 979   if (TMath::Abs(fMagF)>0.005){
 980     //
 981     // use Kalman if mag field
 982     //
 983     ProcessDiff(tp1,tp2, seed,s1,s2);
 984     FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
 985     FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
 986     FillHisto(t1,t2,s1,s2);
 987     ProcessAlign(t1,t2,s1,s2);
 988   }
 989 }
 990
 991 void AliTPCcalibAlign::ProcessAlign(Double_t * t1,
 992                                     Double_t * t2,
 993                                     Int_t s1,Int_t s2){
 994   //
 995   // Do intersector alignment
 996   //
 997   //Process12(t1,t2,GetOrMakeFitter12(s1,s2));
 998   //Process9(t1,t2,GetOrMakeFitter9(s1,s2));
 999   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
1000   ++fPoints[GetIndex(s1,s2)];
1001 }
1002
1003
1004
1005 Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
1006                                        const AliExternalTrackParam &p2) const
1007 {
1008
1009   //
1010   // Accept pair of tracklets?
1011   //
1012   /*
1013   // resolution cuts
1014   TCut cutS0("sqrt(tp2.fC[0]+tp1.fC[0])<0.2");
1015   TCut cutS1("sqrt(tp2.fC[2]+tp1.fC[2])<0.2");
1016   TCut cutS2("sqrt(tp2.fC[5]+tp1.fC[5])<0.01");
1017   TCut cutS3("sqrt(tp2.fC[9]+tp1.fC[9])<0.01");
1018   TCut cutS4("sqrt(tp2.fC[14]+tp1.fC[14])<0.25");
1019   TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
1020   //
1021   // parameters matching cuts
1022   TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1023   TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1024   TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03");
1025   TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03");
1026   TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
1027   TCut cutPP4("abs(tp1.fP[4]-tp2.fP[4])/sqrt(tp2.fC[14]+tp1.fC[14])<3");
1028   TCut cutP=cutP0+cutP1+cutP2+cutP3+cutP4+cutPP4;
1029   */
1030   //
1031   // resolution cuts
1032   Int_t reject=0;
1033   const Double_t *cp1 = p1.GetCovariance();
1034   const Double_t *cp2 = p2.GetCovariance();
1035   if (TMath::Sqrt(cp1[0]+cp2[0])>0.2)  reject|=1;;
1036   if (TMath::Sqrt(cp1[2]+cp2[2])>0.2)  reject|=2;
1037   if (TMath::Sqrt(cp1[5]+cp2[5])>0.01) reject|=4;
1038   if (TMath::Sqrt(cp1[9]+cp2[9])>0.01) reject|=8;
1039   if (TMath::Sqrt(cp1[14]+cp2[14])>0.2) reject|=16;
1040
1041   //parameters difference
1042   const Double_t *tp1 = p1.GetParameter();
1043   const Double_t *tp2 = p2.GetParameter();
1044   if (TMath::Abs(tp1[0]-tp2[0])>0.6) reject|=32;
1045   if (TMath::Abs(tp1[1]-tp2[1])>0.6) reject|=64;
1046   if (TMath::Abs(tp1[2]-tp2[2])>0.03) reject|=128;
1047   if (TMath::Abs(tp1[3]-tp2[3])>0.03) reject|=526;
1048   if (TMath::Abs(tp1[4]-tp2[4])>0.4) reject|=1024;
1049   if (TMath::Abs(tp1[4]-tp2[4])/TMath::Sqrt(cp1[14]+cp2[14])>4) reject|=2048;
1050
1051   //
1052   if (TMath::Abs(tp2[1])>235) reject|=2*4096;
1053
1054   if (fNoField){
1055
1056   }
1057
1058   return reject;
1059 }
1060
1061
1062 Int_t  AliTPCcalibAlign::AcceptTracklet(const Double_t *t1, const Double_t *t2) const
1063 {
1064   //
1065   // accept tracklet  -
1066   //  dist cut + 6 sigma cut
1067   //
1068   Double_t dy     = t2[1]-t1[1];
1069   Double_t dphi   = t2[2]-t1[2];
1070   Double_t dz     = t2[3]-t1[3];
1071   Double_t dtheta = t2[4]-t1[4];
1072   //
1073   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]+0.05*0.05);
1074   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]+0.001*0.001);
1075   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]+0.05*0.05);
1076   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]+0.001*0.001);
1077   //
1078   Int_t reject=0;
1079   if (TMath::Abs(dy)>1.)         reject|=2;
1080   if (TMath::Abs(dphi)>0.1)      reject|=4;
1081   if (TMath::Abs(dz)>1.)         reject|=8;
1082   if (TMath::Abs(dtheta)>0.1)    reject|=16;
1083   //
1084   if (TMath::Abs(dy/sy)>6)         reject|=32;
1085   if (TMath::Abs(dphi/sdydx)>6)    reject|=64;
1086   if (TMath::Abs(dz/sz)>6)         reject|=128;
1087   if (TMath::Abs(dtheta/sdzdx)>6)  reject|=256;
1088   return reject;
1089 }
1090
1091
1092 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
1093                                     const AliExternalTrackParam &t2,
1094                                     const AliTPCseed *seed,
1095                                     Int_t s1,Int_t s2)
1096 {
1097   //
1098   // Process local residuals function
1099   //
1100   TVectorD vecX(160);
1101   TVectorD vecY(160);
1102   TVectorD vecZ(160);
1103   TVectorD vecClY(160);
1104   TVectorD vecClZ(160);
1105   TClonesArray arrCl("AliTPCclusterMI",160);
1106   arrCl.ExpandCreateFast(160);
1107   Int_t count1=0, count2=0;
1108
1109   for (Int_t i=0;i<160;++i) {
1110     AliTPCclusterMI *c=seed->GetClusterPointer(i);
1111     vecX[i]=0;
1112     vecY[i]=0;
1113     vecZ[i]=0;
1114     if (!c) continue;
1115     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
1116     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
1117     vecClY[i] = c->GetY();
1118     vecClZ[i] = c->GetZ();
1119     cl=*c;
1120     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
1121     if (c->GetDetector()==s1) ++count1;
1122     if (c->GetDetector()==s2) ++count2;
1123     Double_t gxyz[3],xyz[3];
1124     t1.GetXYZ(gxyz);
1125     Float_t bz = AliTracker::GetBz(gxyz);
1126     par->GetYAt(c->GetX(), bz, xyz[1]);
1127     par->GetZAt(c->GetX(), bz, xyz[2]);
1128     vecX[i] = c->GetX();
1129     vecY[i]= xyz[1];
1130     vecZ[i]= xyz[2];
1131   }
1132   //
1133   //
1134   if (fStreamLevel>5 && count1>10 && count2>10){
1135     //
1136     // huge output - cluster residuals to be investigated
1137     //
1138     TTreeSRedirector *cstream = GetDebugStreamer();
1139     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
1140     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
1141     /*
1142
1143       Track->Draw("Cl[].fY-vtY.fElements:vtY.fElements-vtX.fElements*tan(pi/18.)>>his(100,-10,0)","Cl.fY!=0&&abs(Cl.fY-vtY.fElements)<1","prof");
1144
1145     */
1146
1147     if (cstream){
1148       (*cstream)<<"Track"<<
1149         "run="<<fRun<<              //  run number
1150         "event="<<fEvent<<          //  event number
1151         "time="<<fTime<<            //  time stamp of event
1152         "trigger="<<fTrigger<<      //  trigger
1153         "triggerClass="<<&fTriggerClass<<      //  trigger
1154         "mag="<<fMagF<<             //  magnetic field
1155         "Cl.="<<&arrCl<<
1156         //"tp0.="<<p0<<
1157         "tp1.="<<p1<<
1158         "tp2.="<<p2<<
1159         "vtX.="<<&vecX<<
1160         "vtY.="<<&vecY<<
1161         "vtZ.="<<&vecZ<<
1162         "vcY.="<<&vecClY<<
1163         "vcZ.="<<&vecClZ<<
1164         "s1="<<s1<<
1165         "s2="<<s2<<
1166         "c1="<<count1<<
1167         "c2="<<count2<<
1168         "\n";
1169     }
1170   }
1171 }
1172
1173
1174
1175
1176 void AliTPCcalibAlign::Process12(const Double_t *t1,
1177                                  const Double_t *t2,
1178                                  TLinearFitter *fitter) const
1179 {
1180   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1181   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1182   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1183   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1184   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1185   //
1186   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
1187   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
1188   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
1189
1190
1191
1192   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1193   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1194
1195   //
1196   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1197   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1198   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1199   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1200
1201   Double_t p[12];
1202   Double_t value;
1203
1204   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1205   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
1206   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1207   for (Int_t i=0; i<12;i++) p[i]=0.;
1208   p[3+0] = x1;          // a10
1209   p[3+1] = y1;          // a11
1210   p[3+2] = z1;          // a12
1211   p[9+1] = 1.;          // a13
1212   p[0+1] = y1*dydx2;    // a01
1213   p[0+2] = z1*dydx2;    // a02
1214   p[9+0] = dydx2;       // a03
1215   value  = y2;
1216   fitter->AddPoint(p,value,sy);
1217
1218   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1219   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
1220   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1221   for (Int_t i=0; i<12;i++) p[i]=0.;
1222   p[6+0] = x1;           // a20
1223   p[6+1] = y1;           // a21
1224   p[6+2] = z1;           // a22
1225   p[9+2] = 1.;           // a23
1226   p[0+1] = y1*dzdx2;     // a01
1227   p[0+2] = z1*dzdx2;     // a02
1228   p[9+0] = dzdx2;        // a03
1229   value  = z2;
1230   fitter->AddPoint(p,value,sz);
1231
1232   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1233   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1234   for (Int_t i=0; i<12;i++) p[i]=0.;
1235   p[3+0] = 1.;           // a10
1236   p[3+1] = dydx1;        // a11
1237   p[3+2] = dzdx1;        // a12
1238   p[0+0] = -dydx2;       // a00
1239   p[0+1] = -dydx1*dydx2; // a01
1240   p[0+2] = -dzdx1*dydx2; // a02
1241   value  = 0.;
1242   fitter->AddPoint(p,value,sdydx);
1243
1244   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1245   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1246   for (Int_t i=0; i<12;i++) p[i]=0.;
1247   p[6+0] = 1;            // a20
1248   p[6+1] = dydx1;        // a21
1249   p[6+2] = dzdx1;        // a22
1250   p[0+0] = -dzdx2;       // a00
1251   p[0+1] = -dydx1*dzdx2; // a01
1252   p[0+2] = -dzdx1*dzdx2; // a02
1253   value  = 0.;
1254   fitter->AddPoint(p,value,sdzdx);
1255 }
1256
1257 void AliTPCcalibAlign::Process9(const Double_t * const t1,
1258                                 const Double_t * const t2,
1259                                 TLinearFitter *fitter) const
1260 {
1261   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1262   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1263   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1264   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1265   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1266   //
1267   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
1268   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
1269   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8]
1270
1271
1272   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1273   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1274   //
1275   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1276   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1277   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1278   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1279
1280   //
1281   Double_t p[12];
1282   Double_t value;
1283
1284   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1285   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
1286   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1287   for (Int_t i=0; i<12;i++) p[i]=0.;
1288   p[2]   += x1;           // a10
1289   //p[]  +=1;             // a11
1290   p[3]   += z1;           // a12
1291   p[7]   += 1;            // a13
1292   p[0]   += y1*dydx2;     // a01
1293   p[1]   += z1*dydx2;     // a02
1294   p[6]   += dydx2;        // a03
1295   value   = y2-y1;        //-a11
1296   fitter->AddPoint(p,value,sy);
1297   //
1298   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1299   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
1300   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1301   for (Int_t i=0; i<12;i++) p[i]=0.;
1302   p[4]   += x1;           // a20
1303   p[5]   += y1;           // a21
1304   //p[]  += z1;           // a22
1305   p[8]   += 1.;           // a23
1306   p[0]   += y1*dzdx2;     // a01
1307   p[1]   += z1*dzdx2;     // a02
1308   p[6]   += dzdx2;        // a03
1309   value  = z2-z1;         //-a22
1310   fitter->AddPoint(p,value,sz);
1311
1312   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1313   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1314   for (Int_t i=0; i<12;i++) p[i]=0.;
1315   p[2]   += 1.;           // a10
1316   //p[]  += dydx1;      // a11
1317   p[3]   += dzdx1;        // a12
1318   //p[]  += -dydx2;       // a00
1319   p[0]   += -dydx1*dydx2; // a01
1320   p[1]   += -dzdx1*dydx2; // a02
1321   value  = -dydx1+dydx2;  // -a11 + a00
1322   fitter->AddPoint(p,value,sdydx);
1323
1324   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1325   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1326   for (Int_t i=0; i<12;i++) p[i]=0.;
1327   p[4]   += 1;            // a20
1328   p[5]   += dydx1;        // a21
1329   //p[]  += dzdx1;        // a22
1330   //p[]  += -dzdx2;       // a00
1331   p[0]   += -dydx1*dzdx2; // a01
1332   p[1]   += -dzdx1*dzdx2; // a02
1333   value  = -dzdx1+dzdx2;  // -a22 + a00
1334   fitter->AddPoint(p,value,sdzdx);
1335 }
1336
1337 void AliTPCcalibAlign::Process6(const Double_t *const t1,
1338                                 const Double_t *const t2,
1339                                 TLinearFitter *fitter) const
1340 {
1341   // x2    =  1  *x1 +-a01*y1 + 0      +a03
1342   // y2    =  a01*x1 + 1  *y1 + 0      +a13
1343   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
1344   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1345   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1346   //
1347   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
1348   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
1349   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5]
1350
1351   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1352   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1353
1354   //
1355   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1356   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1357   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1358   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1359
1360
1361   Double_t p[12];
1362   Double_t value;
1363   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1364   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
1365   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1366   for (Int_t i=0; i<12;i++) p[i]=0.;
1367   p[0]   += x1;           // a10
1368   //p[]  +=1;             // a11
1369   //p[]  += z1;           // a12
1370   p[4]   += 1;            // a13
1371   p[0]   += -y1*dydx2;    // a01
1372   //p[]  += z1*dydx2;     // a02
1373   p[3]   += dydx2;        // a03
1374   value   = y2-y1;        //-a11
1375   fitter->AddPoint(p,value,sy);
1376   //
1377   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1378   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
1379   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1380   for (Int_t i=0; i<12;i++) p[i]=0.;
1381   p[1]   += x1;           // a20
1382   p[2]   += y1;           // a21
1383   //p[]  += z1;           // a22
1384   p[5]   += 1.;           // a23
1385   p[0]   += -y1*dzdx2;    // a01
1386   //p[]   += z1*dzdx2;     // a02
1387   p[3]   += dzdx2;        // a03
1388   value  = z2-z1;         //-a22
1389   fitter->AddPoint(p,value,sz);
1390
1391   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1392   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1393   for (Int_t i=0; i<12;i++) p[i]=0.;
1394   p[0]   += 1.;           // a10
1395   //p[]  += dydx1;      // a11
1396   //p[]   += dzdx1;        // a12
1397   //p[]  += -dydx2;       // a00
1398   //p[0]   +=  dydx1*dydx2; // a01         FIXME- 0912 MI
1399   //p[]   += -dzdx1*dydx2; // a02
1400   value  = -dydx1+dydx2;  // -a11 + a00
1401   fitter->AddPoint(p,value,sdydx);
1402
1403   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1404   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1405   for (Int_t i=0; i<12;i++) p[i]=0.;
1406   p[1]   += 1;            // a20
1407   //  p[2]   += dydx1;        // a21   FIXME- 0912 MI
1408   //p[]  += dzdx1;        // a22
1409   //p[]  += -dzdx2;       // a00
1410   //p[0]   +=  dydx1*dzdx2; // a01     FIXME- 0912 MI
1411   //p[]  += -dzdx1*dzdx2; // a02
1412   value  = -dzdx1+dzdx2;  // -a22 + a00
1413   fitter->AddPoint(p,value,sdzdx);
1414 }
1415
1416
1417
1418
1419 void AliTPCcalibAlign::EvalFitters(Int_t minPoints) {
1420   //
1421   // Analyze function
1422   //
1423   // Perform the fitting using linear fitters
1424   //
1425   TLinearFitter *f;
1426   TFile fff("alignDebug.root","recreate");
1427   for (Int_t s1=0;s1<72;++s1)
1428     for (Int_t s2=0;s2<72;++s2){
1429       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1430         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1431         if (f->Eval()!=0) {
1432           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1433           f->Write(Form("f12_%d_%d",s1,s2));
1434         }else{
1435           f->Write(Form("f12_%d_%d",s1,s2));
1436         }
1437       }
1438       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1439         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1440         if (f->Eval()!=0) {
1441           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1442         }else{
1443           f->Write(Form("f9_%d_%d",s1,s2));
1444         }
1445       }
1446       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1447         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1448         if (f->Eval()!=0) {
1449           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1450         }else{
1451           f->Write(Form("f6_%d_%d",s1,s2));
1452         }
1453       }
1454     }
1455   TMatrixD mat(4,4);
1456   for (Int_t s1=0;s1<72;++s1)
1457     for (Int_t s2=0;s2<72;++s2){
1458       if (GetTransformation12(s1,s2,mat)){
1459         fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2));
1460       }
1461       if (GetTransformation9(s1,s2,mat)){
1462         fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2));
1463       }
1464       if (GetTransformation6(s1,s2,mat)){
1465         fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2));
1466       }
1467     }
1468   //this->Write("align");
1469
1470 }
1471
1472 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
1473   //
1474   // get or make fitter - general linear transformation
1475   //
1476   static Int_t counter12=0;
1477   static TF1 f12("f12","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]++x[9]++x[10]++x[11]");
1478   TLinearFitter * fitter = GetFitter12(s1,s2);
1479   if (fitter) return fitter;
1480   //  fitter =new TLinearFitter(12,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]++x[9]++x[10]++x[11]");
1481   fitter =new TLinearFitter(&f12,"");
1482   fitter->StoreData(kFALSE);
1483   fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
1484   counter12++;
1485   //  if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
1486   return fitter;
1487 }
1488
1489 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
1490   //
1491   //get or make fitter - general linear transformation - no scaling
1492   //
1493   static Int_t counter9=0;
1494   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
1495   TLinearFitter * fitter = GetFitter9(s1,s2);
1496   if (fitter) return fitter;
1497   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
1498   fitter =new TLinearFitter(&f9,"");
1499   fitter->StoreData(kFALSE);
1500   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
1501   counter9++;
1502   //  if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
1503   return fitter;
1504 }
1505
1506 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
1507   //
1508   // get or make fitter  - 6 paramater linear tranformation
1509   //                     - no scaling
1510   //                     - rotation x-y
1511   //                     - tilting x-z, y-z
1512   static Int_t counter6=0;
1513   //static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");  // change in the back compatibility in root -
1514   TLinearFitter * fitter = GetFitter6(s1,s2);
1515   if (fitter) return fitter;
1516   fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1517   //fitter=new TLinearFitter(&f6,"");
1518   fitter->StoreData(kFALSE);
1519   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
1520   counter6++;
1521   //  if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
1522   return fitter;
1523   /*
1524   //    as root changed interface at some moment I put here also the test - which should be unit test in future:
1525   TLinearFitter fitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1526   for (Int_t i=0; i<10000; i++) {
1527     x[0]=i;
1528     x[1]=i*i;
1529     x[2]=i*i*i;
1530     x[3]=TMath::Power(i,1/2.);
1531     x[4]=TMath::Power(i,1/3.);
1532     x[5]=TMath::Power(i,1/4.);
1533     fitter.AddPoint(x,i*i*i);
1534   }
1535   fitter.Eval();
1536   fitter.GetParameters(vec);
1537   vec.Print();
1538   => result of the fit should = for all vector elements except of element 2; - it was the case with root v5-34-09
1539    */
1540 }
1541
1542
1543
1544
1545
1546 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
1547   //
1548   // GetTransformation matrix - 12 paramaters - generael linear transformation
1549   //
1550   if (!GetFitter12(s1,s2))
1551     return false;
1552   else {
1553     TVectorD p(12);
1554     GetFitter12(s1,s2)->GetParameters(p);
1555     a.ResizeTo(4,4);
1556     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
1557     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
1558     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
1559     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1560     return true;
1561   }
1562 }
1563
1564 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
1565   //
1566   // GetTransformation matrix - 9 paramaters - general linear transformation
1567   //                            No scaling
1568   //
1569   if (!GetFitter9(s1,s2))
1570     return false;
1571   else {
1572     TVectorD p(9);
1573     GetFitter9(s1,s2)->GetParameters(p);
1574     a.ResizeTo(4,4);
1575     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
1576     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
1577     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
1578     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1579     return true;
1580   }
1581 }
1582
1583 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
1584   //
1585   // GetTransformation matrix - 6  paramaters
1586   //                            3  translation
1587   //                            1  rotation -x-y
1588   //                            2  tilting x-z y-z
1589   if (!GetFitter6(s1,s2))
1590     return false;
1591   else {
1592     TVectorD p(6);
1593     GetFitter6(s1,s2)->GetParameters(p);
1594     a.ResizeTo(4,4);
1595     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
1596     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
1597     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
1598     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1599     return true;
1600   }
1601 }
1602
1603 void AliTPCcalibAlign::MakeResidualHistos(){
1604   //
1605   // Make cluster residual histograms
1606   //
1607   Double_t xminTrack[9], xmaxTrack[9];
1608   Int_t    binsTrack[9];
1609   TString  axisName[9],axisTitle[9];
1610   //
1611   // 0 - delta   of interest
1612   // 1 - global  phi in sector number  as float
1613   // 2 - local   x
1614   // 3 - local   ky
1615   // 4 - local   kz
1616   //
1617   axisName[0]="delta";   axisTitle[0]="#Delta (cm)";
1618   if (TMath::Abs(AliTracker::GetBz())<0.01){
1619     binsTrack[0]=60;       xminTrack[0]=-1.2;        xmaxTrack[0]=1.2;
1620   }else{
1621     binsTrack[0]=60;       xminTrack[0]=-0.6;        xmaxTrack[0]=0.6;
1622   }
1623   //
1624   axisName[1]="sector";   axisTitle[1]="Sector Number";
1625   binsTrack[1]=180;       xminTrack[1]=0;        xmaxTrack[1]=18;
1626   //
1627   axisName[2]="R";   axisTitle[2]="r (cm)";
1628   binsTrack[2]=53;       xminTrack[2]=85.;        xmaxTrack[2]=245.;
1629   //
1630   //
1631   axisName[3]="kZ";      axisTitle[3]="dz/dx";
1632   binsTrack[3]=36;       xminTrack[3]=-1.8;        xmaxTrack[3]=1.8;
1633   //
1634   fClusterDelta[0] = new THnF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1635   fClusterDelta[1] = new THnF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1636   //
1637   //
1638   //
1639   for (Int_t ivar=0;ivar<2;ivar++){
1640     for (Int_t ivar2=0;ivar2<4;ivar2++){
1641       fClusterDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1642       fClusterDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
1643     }
1644   }
1645
1646 }
1647
1648
1649 void AliTPCcalibAlign::MakeResidualHistosTracklet(){
1650   //
1651   // Make tracklet residual histograms
1652   //
1653   Double_t xminTrack[9], xmaxTrack[9];
1654   Int_t    binsTrack[9];
1655   TString  axisName[9],axisTitle[9];
1656   //
1657   // 0 - delta   of interest
1658   // 1 - global  phi in sector number  as float
1659   // 2 - local   x
1660   // 3 - local   ky
1661   // 4 - local   kz
1662   // 5 - sector  1
1663   // 6 - sector  0
1664   // 7 - z position  0
1665
1666   axisName[0]="delta";   axisTitle[0]="#Delta (cm)";
1667   binsTrack[0]=60;       xminTrack[0]=-0.5;        xmaxTrack[0]=0.5;
1668   //
1669   axisName[1]="phi";   axisTitle[1]="#phi";
1670   binsTrack[1]=90;       xminTrack[1]=-TMath::Pi();        xmaxTrack[1]=TMath::Pi();
1671   //
1672   axisName[2]="localX";   axisTitle[2]="x (cm)";
1673   binsTrack[2]=10;       xminTrack[2]=120.;        xmaxTrack[2]=200.;
1674   //
1675   axisName[3]="kY";      axisTitle[3]="dy/dx";
1676   binsTrack[3]=10;       xminTrack[3]=-0.5;        xmaxTrack[3]=0.5;
1677   //
1678   axisName[4]="kZ";      axisTitle[4]="dz/dx";
1679   binsTrack[4]=11;       xminTrack[4]=-1.1;        xmaxTrack[4]=1.1;
1680   //
1681   axisName[5]="is1";      axisTitle[5]="is1";
1682   binsTrack[5]=72;       xminTrack[5]=0;        xmaxTrack[5]=72;
1683   //
1684   axisName[6]="is0";      axisTitle[6]="is0";
1685   binsTrack[6]=72;       xminTrack[6]=0;        xmaxTrack[6]=72;
1686   //
1687   axisName[7]="z";        axisTitle[7]="z(cm)";
1688   binsTrack[7]=12;        xminTrack[7]=-240;        xmaxTrack[7]=240;
1689   //
1690   axisName[8]="IsPrimary";        axisTitle[8]="Is Primary";
1691   binsTrack[8]=2;        xminTrack[8]=-0.1;        xmaxTrack[8]=1.1;
1692
1693   //
1694   xminTrack[0]=-0.25;        xmaxTrack[0]=0.25;
1695   fTrackletDelta[0] = new THnSparseF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
1696   xminTrack[0]=-0.5;        xmaxTrack[0]=0.5;
1697   fTrackletDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
1698   xminTrack[0]=-0.005;        xmaxTrack[0]=0.005;
1699   fTrackletDelta[2] = new THnSparseF("#Delta_{kY}","#Delta_{kY}", 9, binsTrack,xminTrack, xmaxTrack);
1700   xminTrack[0]=-0.008;        xmaxTrack[0]=0.008;
1701   fTrackletDelta[3] = new THnSparseF("#Delta_{kZ}","#Delta_{kZ}", 9, binsTrack,xminTrack, xmaxTrack);
1702   //
1703   //
1704   //
1705   for (Int_t ivar=0;ivar<4;ivar++){
1706     for (Int_t ivar2=0;ivar2<9;ivar2++){
1707       fTrackletDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1708       fTrackletDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
1709     }
1710   }
1711
1712 }
1713
1714
1715
1716 void AliTPCcalibAlign::FillHisto(const Double_t *t1,
1717                                  const Double_t *t2,
1718                                  Int_t s1,Int_t s2) {
1719   //
1720   // Fill residual histograms
1721   // Track2-Track1
1722   // Innner-Outer
1723   // Left right - x-y
1724   // A-C side
1725   if (1)  {
1726     Double_t dy     = t2[1]-t1[1];
1727     Double_t dphi   = t2[2]-t1[2];
1728     Double_t dz     = t2[3]-t1[3];
1729     Double_t dtheta = t2[4]-t1[4];
1730     Double_t zmean = (t2[3]+t1[3])*0.5;
1731     //
1732     GetHisto(kPhi,s1,s2,kTRUE)->Fill(dphi);
1733     GetHisto(kTheta,s1,s2,kTRUE)->Fill(dtheta);
1734     GetHisto(kY,s1,s2,kTRUE)->Fill(dy);
1735     GetHisto(kZ,s1,s2,kTRUE)->Fill(dz);
1736     //
1737     GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(zmean,dphi);
1738     GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(zmean,dtheta);
1739     GetHisto(kYz,s1,s2,kTRUE)->Fill(zmean,dy);
1740     GetHisto(kZz,s1,s2,kTRUE)->Fill(zmean,dz);
1741     //
1742     GetHisto(kYPhi,s1,s2,kTRUE)->Fill(t2[2],dy);
1743     GetHisto(kZTheta,s1,s2,kTRUE)->Fill(t2[4],dz);
1744   }
1745 }
1746
1747
1748 void AliTPCcalibAlign::FillHisto(AliExternalTrackParam *tp1,
1749                                  AliExternalTrackParam *tp2,
1750                                  Int_t s1,Int_t s2) {
1751   //
1752   // Fill residual histograms
1753   // Track2-Track1
1754   if (s2<s1) return;//
1755   const Double_t kEpsilon=0.001;
1756   Double_t x[9]={0,0,0,0,0,0,0,0,0};
1757   AliExternalTrackParam p1(*tp1);
1758   AliExternalTrackParam p2(*tp2);
1759   if (s1%18==s2%18) {
1760     // inner outer - match at the IROC-OROC boundary
1761     if (!p1.PropagateTo(fXIO, AliTrackerBase::GetBz())) return;
1762   }
1763   if (!p2.Rotate(p1.GetAlpha())) return;
1764   if (!p2.PropagateTo(p1.GetX(),AliTrackerBase::GetBz())) return;
1765   if (TMath::Abs(p1.GetX()-p2.GetX())>kEpsilon) return;
1766   Double_t xyz[3];
1767   p1.GetXYZ(xyz);
1768   x[1]=TMath::ATan2(xyz[1],xyz[0]);
1769   x[2]=p1.GetX();
1770   x[3]=0.5*(p1.GetSnp()+p2.GetSnp());  // mean snp
1771   x[4]=0.5*(p1.GetTgl()+p2.GetTgl());  // mean tgl
1772   x[5]=s2;
1773   x[6]=s1;
1774   x[7]=0.5*(p1.GetZ()+p2.GetZ());
1775   // is primary ?
1776   Int_t  isPrimary = (TMath::Abs(p1.GetTgl()-p1.GetZ()/p1.GetX())<0.1) ? 1:0;
1777   x[8]= isPrimary;
1778   //
1779   x[0]=p2.GetY()-p1.GetY();
1780   fTrackletDelta[0]->Fill(x);
1781   x[0]=p2.GetZ()-p1.GetZ();
1782   fTrackletDelta[1]->Fill(x);
1783   x[0]=p2.GetSnp()-p1.GetSnp();
1784   fTrackletDelta[2]->Fill(x);
1785   x[0]=p2.GetTgl()-p1.GetTgl();
1786   fTrackletDelta[3]->Fill(x);
1787   TTreeSRedirector *cstream = GetDebugStreamer();
1788   if (cstream){
1789     (*cstream)<<"trackletMatch"<<
1790       "tp1.="<<tp1<<   // input tracklet
1791       "tp2.="<<tp2<<
1792       "p1.="<<&p1<<    // tracklet in the ref frame
1793       "p2.="<<&p2<<
1794       "s1="<<s1<<
1795       "s2="<<s2<<
1796       "\n";
1797   }
1798
1799 }
1800
1801
1802
1803 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
1804 {
1805   //
1806   // return specified residual histogram - it is only QA
1807   // if force specified the histogram and given histogram is not existing
1808   //  new histogram is created
1809   //
1810   if (GetIndex(s1,s2)>=72*72) return 0;
1811   TObjArray *histoArray=0;
1812   switch (type) {
1813   case kY:
1814     histoArray = &fDyHistArray; break;
1815   case kZ:
1816     histoArray = &fDzHistArray; break;
1817   case kPhi:
1818     histoArray = &fDphiHistArray; break;
1819   case kTheta:
1820     histoArray = &fDthetaHistArray; break;
1821   case kYPhi:
1822     histoArray = &fDyPhiHistArray; break;
1823   case kZTheta:
1824     histoArray = &fDzThetaHistArray; break;
1825   case kYz:
1826     histoArray = &fDyZHistArray; break;
1827   case kZz:
1828     histoArray = &fDzZHistArray; break;
1829   case kPhiZ:
1830     histoArray = &fDphiZHistArray; break;
1831   case kThetaZ:
1832     histoArray = &fDthetaZHistArray; break;
1833   }
1834   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
1835   if (histo) return histo;
1836   if (force==kFALSE) return 0;
1837   //
1838   stringstream name;
1839   stringstream title;
1840   switch (type) {
1841   case kY:
1842     name<<"hist_y_"<<s1<<"_"<<s2;
1843     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1844     histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.5,0.5); // +/- 5 mm
1845     break;
1846   case kZ:
1847     name<<"hist_z_"<<s1<<"_"<<s2;
1848     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1849     histo = new TH1D(name.str().c_str(),title.str().c_str(),100,-0.3,0.3); // +/- 3 mm
1850     break;
1851   case kPhi:
1852     name<<"hist_phi_"<<s1<<"_"<<s2;
1853     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1854     histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.01,0.01); // +/- 10 mrad
1855     break;
1856   case kTheta:
1857     name<<"hist_theta_"<<s1<<"_"<<s2;
1858     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1859     histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.01,0.01); // +/- 10 mrad
1860     break;
1861     //
1862     //
1863   case kYPhi:
1864     name<<"hist_yphi_"<<s1<<"_"<<s2;
1865     title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1866     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,100,-0.5,0.5); // +/- 5 mm
1867     break;
1868   case kZTheta:
1869     name<<"hist_ztheta_"<<s1<<"_"<<s2;
1870     title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1871     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,100,-0.3,0.3); // +/- 3 mm
1872     break;
1873     //
1874     //
1875     //
1876   case kYz:
1877     name<<"hist_yz_"<<s1<<"_"<<s2;
1878     title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1879     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.5,0.5); // +/- 5 mm
1880     break;
1881   case kZz:
1882     name<<"hist_zz_"<<s1<<"_"<<s2;
1883     title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1884     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.3,0.3); // +/- 3 mm
1885     break;
1886   case kPhiZ:
1887     name<<"hist_phiz_"<<s1<<"_"<<s2;
1888     title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1889     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.01,0.01); // +/- 10 mrad
1890     break;
1891   case kThetaZ:
1892     name<<"hist_thetaz_"<<s1<<"_"<<s2;
1893     title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1894     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.01,0.01); // +/- 10 mrad
1895     break;
1896
1897
1898   }
1899   histo->SetDirectory(0);
1900   histoArray->AddAt(histo,GetIndex(s1,s2));
1901   return histo;
1902 }
1903
1904 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
1905                                            Int_t i0, Int_t i1, FitType type)
1906 {
1907   //
1908   //
1909   //
1910   TMatrixD mat;
1911   //TObjArray *fitArray=0;
1912   Double_t xsec[1000];
1913   Double_t ysec[1000];
1914   Int_t npoints=0;
1915   for (Int_t isec = sec0; isec<=sec1; isec++){
1916     Int_t isec2 = (isec+dsec)%72;
1917     switch (type) {
1918     case k6:
1919       GetTransformation6(isec,isec2,mat);break;
1920     case k9:
1921       GetTransformation9(isec,isec2,mat);break;
1922     case k12:
1923       GetTransformation12(isec,isec2,mat);break;
1924     }
1925     xsec[npoints]=isec;
1926     ysec[npoints]=mat(i0,i1);
1927     ++npoints;
1928   }
1929   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1930   Char_t name[1000];
1931   snprintf(name,100,"Mat[%d,%d]  Type=%d",i0,i1,type);
1932   gr->SetName(name);
1933   return gr;
1934 }
1935
1936 void  AliTPCcalibAlign::MakeTree(const char *fname, Int_t minPoints){
1937   //
1938   // make tree with alignment cosntant  -
1939   // For  QA visualization
1940   //
1941   /*
1942     TFile fcalib("CalibObjects.root");
1943     TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
1944     AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
1945     align->EvalFitters();
1946     align->MakeTree("alignTree.root");
1947     TFile falignTree("alignTree.root");
1948     TTree * treeAlign = (TTree*)falignTree.Get("Align");
1949    */
1950   TTreeSRedirector cstream(fname);
1951   for (Int_t s1=0;s1<72;++s1)
1952     for (Int_t s2=0;s2<72;++s2){
1953       TMatrixD m6;
1954       TMatrixD m6FX;
1955       TMatrixD m9;
1956       TMatrixD m12;
1957       TVectorD param6Diff;  // align parameters diff
1958       TVectorD param6s1(6);    // align parameters sector1
1959       TVectorD param6s2(6);    // align parameters sector2
1960
1961       //
1962       //
1963       if (fSectorParamA){
1964         TMatrixD * kpar = fSectorParamA;
1965         TMatrixD * kcov = fSectorCovarA;
1966         if (s1%36>=18){
1967           kpar = fSectorParamC;
1968           kcov = fSectorCovarC;
1969         }
1970         for (Int_t ipar=0;ipar<6;ipar++){
1971           Int_t isec1 = s1%18;
1972           Int_t isec2 = s2%18;
1973           if (s1>35) isec1+=18;
1974           if (s2>35) isec2+=18;
1975           param6s1(ipar)=(*kpar)(6*isec1+ipar,0);
1976           param6s2(ipar)=(*kpar)(6*isec2+ipar,0);
1977         }
1978       }
1979
1980       Double_t dy=0, dz=0, dphi=0,dtheta=0;
1981       Double_t sy=0, sz=0, sphi=0,stheta=0;
1982       Double_t ny=0, nz=0, nphi=0,ntheta=0;
1983       Double_t chi2v12=0, chi2v9=0, chi2v6=0;
1984       //      Int_t npoints=0;
1985       // TLinearFitter * fitter = 0;
1986       if (fPoints[GetIndex(s1,s2)]>minPoints){
1987         //
1988         //
1989         //
1990 //      fitter = GetFitter12(s1,s2);
1991 //      npoints = fitter->GetNpoints();
1992 //      chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1993
1994 //      //
1995 //      fitter = GetFitter9(s1,s2);
1996 //      npoints = fitter->GetNpoints();
1997 //      chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1998 //      //
1999 //      fitter = GetFitter6(s1,s2);
2000 //      npoints = fitter->GetNpoints();
2001 //      chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
2002 //      fitter->GetParameters(param6Diff);
2003 //      //
2004 //      GetTransformation6(s1,s2,m6);
2005 //      GetTransformation9(s1,s2,m9);
2006 //      GetTransformation12(s1,s2,m12);
2007 //      //
2008 //      fitter = GetFitter6(s1,s2);
2009 //      //fitter->FixParameter(3,0);
2010 //      //fitter->Eval();
2011 //      GetTransformation6(s1,s2,m6FX);
2012         //
2013         TH1 * his=0;
2014         his = GetHisto(kY,s1,s2);
2015         if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
2016         his = GetHisto(kZ,s1,s2);
2017         if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
2018         his = GetHisto(kPhi,s1,s2);
2019         if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
2020         his = GetHisto(kTheta,s1,s2);
2021         if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
2022         //
2023
2024       }
2025       AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2026       if (!magF) AliError("Magneticd field - not initialized");
2027       Double_t bz = magF->SolenoidField()/10.; //field in T
2028
2029       cstream<<"Align"<<
2030         "run="<<fRun<<  // run
2031         "bz="<<bz<<
2032         "s1="<<s1<<     // reference sector
2033         "s2="<<s2<<     // sector to align
2034         "m6FX.="<<&m6FX<<   // tranformation matrix
2035         "m6.="<<&m6<<   // tranformation matrix
2036         "m9.="<<&m9<<   //
2037         "m12.="<<&m12<<
2038         "chi2v12="<<chi2v12<<
2039         "chi2v9="<<chi2v9<<
2040         "chi2v6="<<chi2v6<<
2041         //
2042         "p6.="<<&param6Diff<<
2043         "p6s1.="<<&param6s1<<
2044         "p6s2.="<<&param6s2<<
2045         //               histograms mean RMS and entries
2046         "dy="<<dy<<
2047         "sy="<<sy<<
2048         "ny="<<ny<<
2049         "dz="<<dz<<
2050         "sz="<<sz<<
2051         "nz="<<nz<<
2052         "dphi="<<dphi<<
2053         "sphi="<<sphi<<
2054         "nphi="<<nphi<<
2055         "dtheta="<<dtheta<<
2056         "stheta="<<stheta<<
2057         "ntheta="<<ntheta<<
2058         "\n";
2059     }
2060
2061 }
2062
2063
2064 //_____________________________________________________________________
2065 Long64_t AliTPCcalibAlign::Merge(TCollection* const list) {
2066   //
2067   // merge function
2068   //
2069   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
2070   if (!list)
2071     return 0;
2072   if (list->IsEmpty())
2073     return 1;
2074
2075   TIterator* iter = list->MakeIterator();
2076   TObject* obj = 0;
2077   iter->Reset();
2078   Int_t count=0;
2079   //
2080   TString str1(GetName());
2081   while((obj = iter->Next()) != 0)
2082     {
2083       AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
2084       if (entry == 0) continue;
2085       if (str1.CompareTo(entry->GetName())!=0) continue;
2086       Add(entry);
2087       count++;
2088     }
2089   return count;
2090 }
2091
2092
2093 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
2094   //
2095   // Add entry - used for merging of compoents
2096   //
2097   for (Int_t i=0; i<72;i++){
2098     for (Int_t j=0; j<72;j++){
2099       if (align->fPoints[GetIndex(i,j)]<1) continue;
2100       fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
2101       //
2102       //
2103       //
2104       for (Int_t itype=0; itype<10; itype++){
2105         TH1 * his0=0, *his1=0;
2106         his0 = GetHisto((HistoType)itype,i,j);
2107         his1 = align->GetHisto((HistoType)itype,i,j);
2108         if (his1){
2109           if (his0) his0->Add(his1);
2110           else {
2111             his0 = GetHisto((HistoType)itype,i,j,kTRUE);
2112             his0->Add(his1);
2113           }
2114         }
2115       }
2116     }
2117   }
2118   TLinearFitter *f0=0;
2119   TLinearFitter *f1=0;
2120   for (Int_t i=0; i<72;i++){
2121     for (Int_t j=0; j<72;j++){
2122       if (align->fPoints[GetIndex(i,j)]<1) continue;
2123       //
2124       //
2125       // fitter12
2126       f0 =  GetFitter12(i,j);
2127       f1 =  align->GetFitter12(i,j);
2128       if (f1){
2129         if (f0) f0->Add(f1);
2130         else {
2131           fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
2132         }
2133       }
2134       //
2135       // fitter9
2136       f0 =  GetFitter9(i,j);
2137       f1 =  align->GetFitter9(i,j);
2138       if (f1){
2139         if (f0) f0->Add(f1);
2140         else {
2141           fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
2142         }
2143       }
2144       f0 =  GetFitter6(i,j);
2145       f1 =  align->GetFitter6(i,j);
2146       if (f1){
2147         if (f0) f0->Add(f1);
2148         else {
2149           fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
2150         }
2151       }
2152     }
2153   }
2154   //
2155   // Add Kalman filter
2156   //
2157   for (Int_t i=0;i<36;i++){
2158     TMatrixD *par0 = (TMatrixD*)fArraySectorIntParam.At(i);
2159     if (!par0){
2160       MakeSectorKalman();
2161       par0 = (TMatrixD*)fArraySectorIntParam.At(i);
2162     }
2163     TMatrixD *par1 = (TMatrixD*)align->fArraySectorIntParam.At(i);
2164     if (!par1) continue;
2165     //
2166     TMatrixD *cov0 = (TMatrixD*)fArraySectorIntCovar.At(i);
2167     TMatrixD *cov1 = (TMatrixD*)align->fArraySectorIntCovar.At(i);
2168     UpdateSectorKalman(*par0,*cov0,*par1,*cov1);
2169   }
2170   if (!fSectorParamA){
2171     MakeKalman();
2172   }
2173   if (align->fSectorParamA){
2174     UpdateKalman(*fSectorParamA,*fSectorCovarA,*align->fSectorParamA,*align->fSectorCovarA);
2175     UpdateKalman(*fSectorParamC,*fSectorCovarC,*align->fSectorParamC,*align->fSectorCovarC);
2176   }
2177   if (!fClusterDelta[0]) MakeResidualHistos();
2178
2179   for (Int_t i=0; i<2; i++){
2180     if (align->fClusterDelta[i]){
2181       fClusterDelta[i]->Add(align->fClusterDelta[i]);
2182     }
2183   }
2184
2185
2186   for (Int_t i=0; i<4; i++){
2187     if (!fTrackletDelta[i] && align->fTrackletDelta[i]) {
2188       fTrackletDelta[i]= (THnSparse*)(align->fTrackletDelta[i]->Clone());
2189       continue;
2190     }
2191     if (align->fTrackletDelta[i]) {
2192       if (fTrackletDelta[i]->GetEntries()<fgkMergeEntriesCut){
2193         fTrackletDelta[i]->Add(align->fTrackletDelta[i]);
2194       }
2195     }
2196   }
2197
2198 }
2199
2200 Double_t AliTPCcalibAlign::Correct(Int_t type, Int_t value, Int_t s1, Int_t s2, Double_t x1, Double_t y1, Double_t z1, Double_t dydx1,Double_t dzdx1){
2201   //
2202   // GetTransformed value
2203   //
2204   //
2205   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
2206   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
2207   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
2208   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2209   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2210
2211
2212   const TMatrixD * mat = GetTransformation(s1,s2,type);
2213   if (!mat) {
2214     if (value==0) return x1;
2215     if (value==1) return y1;
2216     if (value==2) return z1;
2217     if (value==3) return dydx1;
2218     if (value==4) return dzdx1;
2219     //
2220     if (value==5) return dydx1;
2221     if (value==6) return dzdx1;
2222     return 0;
2223   }
2224   Double_t valT=0;
2225
2226   if (value==0){
2227     valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
2228   }
2229
2230   if (value==1){
2231     valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
2232   }
2233   if (value==2){
2234     valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
2235   }
2236   if (value==3){
2237     //    dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2238     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
2239     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
2240   }
2241
2242   if (value==4){
2243     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2244     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
2245     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
2246   }
2247   //
2248   if (value==5){
2249     // onlys shift in angle
2250     //    dydx2 =  (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2251     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1;
2252   }
2253
2254   if (value==6){
2255     // only shift in angle
2256     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2257     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1;
2258   }
2259   //
2260   return valT;
2261 }
2262
2263
2264 void  AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){
2265   //
2266   // Update track parameters t1
2267   //
2268   TMatrixD vecXk(5,1);    // X vector
2269   TMatrixD covXk(5,5);    // X covariance
2270   TMatrixD matHk(1,5);    // vector to mesurement
2271   TMatrixD measR(1,1);    // measurement error
2272   //TMatrixD matQk(5,5);    // prediction noise vector
2273   TMatrixD vecZk(1,1);    // measurement
2274   //
2275   TMatrixD vecYk(1,1);    // Innovation or measurement residual
2276   TMatrixD matHkT(5,1);
2277   TMatrixD matSk(1,1);    // Innovation (or residual) covariance
2278   TMatrixD matKk(5,1);    // Optimal Kalman gain
2279   TMatrixD mat1(5,5);     // update covariance matrix
2280   TMatrixD covXk2(5,5);   //
2281   TMatrixD covOut(5,5);
2282   //
2283   Double_t *param1=(Double_t*) track1.GetParameter();
2284   Double_t *covar1=(Double_t*) track1.GetCovariance();
2285
2286   //
2287   // copy data to the matrix
2288   for (Int_t ipar=0; ipar<5; ipar++){
2289     vecXk(ipar,0) = param1[ipar];
2290     for (Int_t jpar=0; jpar<5; jpar++){
2291       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
2292     }
2293   }
2294   //
2295   //
2296   //
2297   vecZk(0,0) = track2.GetParameter()[4];   // 1/pt measurement from track 2
2298   measR(0,0) = track2.GetCovariance()[14];  // 1/pt measurement error
2299   if (noField) {
2300     measR(0,0)*=0.000000001;
2301     vecZk(0,0)=0.;
2302   }
2303   //
2304   matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0;
2305   matHk(0,3)= 0;    matHk(0,4)= 1;           // vector to measurement
2306   //
2307   //
2308   //
2309   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
2310   matHkT=matHk.T(); matHk.T();
2311   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
2312   matSk.Invert();
2313   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
2314   vecXk += matKk*vecYk;                      //  updated vector
2315   mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
2316   covXk2 = (mat1-(matKk*matHk));
2317   covOut =  covXk2*covXk;
2318   //
2319   //
2320   //
2321   // copy from matrix to parameters
2322   if (0) {
2323     covOut.Print();
2324     vecXk.Print();
2325     covXk.Print();
2326     track1.Print();
2327     track2.Print();
2328   }
2329
2330   for (Int_t ipar=0; ipar<5; ipar++){
2331     param1[ipar]= vecXk(ipar,0) ;
2332     for (Int_t jpar=0; jpar<5; jpar++){
2333       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
2334     }
2335   }
2336
2337 }
2338
2339 void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){
2340   //
2341   //  Global Align -combine the partial alignment of pair of sectors
2342   //  minPoints - minimal number of points - don't use sector alignment wit smaller number
2343   //  sysError  - error added to the alignemnt error
2344   //
2345   AliTPCcalibAlign * align = this;
2346   TMatrixD * arrayAlign[72];
2347   TMatrixD * arrayAlignDiff[72];
2348   //
2349   for (Int_t i=0;i<72; i++) {
2350     TMatrixD * mat = new TMatrixD(4,4);
2351     mat->UnitMatrix();
2352     arrayAlign[i]=mat;
2353     arrayAlignDiff[i]=(TMatrixD*)(mat->Clone());
2354   }
2355
2356   TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root");
2357   for (Int_t iter=0; iter<niter;iter++){
2358     printf("Iter=\t%d\n",iter);
2359     for (Int_t is0=0;is0<72; is0++) {
2360       //
2361       //TMatrixD  *mati0 = arrayAlign[is0];
2362       TMatrixD matDiff(4,4);
2363       Double_t sumw=0;
2364       for (Int_t is1=0;is1<72; is1++) {
2365         Bool_t invers=kFALSE;
2366         Int_t npoints=0;
2367         TMatrixD covar;
2368         TVectorD errors;
2369         const TMatrixD *mat = align->GetTransformation(is0,is1,0);
2370         if (mat){
2371           npoints = align->GetFitter6(is0,is1)->GetNpoints();
2372           if (npoints>minPoints){
2373             align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar);
2374             align->GetFitter6(is0,is1)->GetErrors(errors);
2375           }
2376         }
2377         else{
2378           invers=kTRUE;
2379           mat = align->GetTransformation(is1,is0,0);
2380           if (mat) {
2381             npoints = align->GetFitter6(is1,is0)->GetNpoints();
2382             if (npoints>minPoints){
2383               align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar);
2384               align->GetFitter6(is1,is0)->GetErrors(errors);
2385             }
2386           }
2387         }
2388         if (!mat) continue;
2389         if (npoints<minPoints) continue;
2390         //
2391         Double_t weight=1;
2392         if (is1/36>is0/36) weight*=2./3.; //IROC-OROC
2393         if (is1/36<is0/36) weight*=1./3.; //OROC-IROC
2394         if (is1/36==is0/36) weight*=1/3.; //OROC-OROC
2395         if (is1%36!=is0%36) weight*=1/2.; //Not up-down
2396         weight/=(errors[4]*errors[4]+sysError*sysError);  // wieghting with error in Y
2397         //
2398         //
2399         TMatrixD matT = *mat;
2400         if (invers)  matT.Invert();
2401         TMatrixD diffMat= (*(arrayAlign[is1]))*matT;
2402         diffMat-=(*arrayAlign[is0]);
2403         matDiff+=weight*diffMat;
2404         sumw+=weight;
2405
2406         (*cstream)<<"LAlign"<<
2407           "iter="<<iter<<
2408           "s0="<<is0<<
2409           "s1="<<is1<<
2410           "npoints="<<npoints<<
2411           "m60.="<<arrayAlign[is0]<<
2412           "m61.="<<arrayAlign[is1]<<
2413           "m01.="<<&matT<<
2414           "diff.="<<&diffMat<<
2415           "cov.="<<&covar<<
2416           "err.="<<&errors<<
2417           "\n";
2418       }
2419       if (sumw>0){
2420         matDiff*=1/sumw;
2421         matDiff(0,0)=0;
2422         matDiff(1,1)=0;
2423         matDiff(1,1)=0;
2424         matDiff(1,1)=0;
2425         (*arrayAlignDiff[is0]) = matDiff;
2426       }
2427     }
2428     for (Int_t is0=0;is0<72; is0++) {
2429       if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]);
2430       if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]);
2431        //
2432       (*cstream)<<"GAlign"<<
2433         "iter="<<iter<<
2434         "s0="<<is0<<
2435         "m6.="<<arrayAlign[is0]<<
2436         "\n";
2437     }
2438   }
2439
2440   delete cstream;
2441   for (Int_t isec=0;isec<72;isec++){
2442     fCombinedMatrixArray6.AddAt(arrayAlign[isec],isec);
2443     delete arrayAlignDiff[isec];
2444   }
2445 }
2446
2447
2448  Int_t  AliTPCcalibAlign::RefitLinear(const AliTPCseed * track, Int_t isec, Double_t *fitParam, Int_t refSector,  TMatrixD &tparam, TMatrixD&tcovar, Double_t xRef, Bool_t both){
2449   //
2450   // Refit tracklet linearly using clusters at  given sector isec
2451   // Clusters are rotated to the  reference frame of sector refSector
2452   //
2453   // fit parameters and errors retruning in the fitParam
2454   //
2455   // seed      - acces to the original clusters
2456   // isec      - sector to be refited
2457   // fitParam  -
2458   //           0  lx
2459   //           1  ly
2460   //           2  dy/dz
2461   //           3  lz
2462   //           4  dz/dx
2463   //           5  sx
2464   //           6  sy
2465   //           7  sdydx
2466   //           8  sz
2467   //           9  sdzdx
2468   // ref sector is the sector defining ref frame - rotation
2469   // return value - number of used clusters
2470
2471   const Int_t      kMinClusterF=15;
2472   const  Int_t     kdrow1 =10;        // rows to skip at the end
2473   const  Int_t     kdrow0 =3;        // rows to skip at beginning
2474   const  Float_t   kedgeyIn=2.5;
2475   const  Float_t   kedgeyOut=4.0;
2476   const  Float_t   kMaxDist=5;       // max distance -in sigma
2477   const  Float_t   kMaxCorrY=0.05;    // max correction
2478   //
2479   Double_t dalpha = 0;
2480   if ((refSector%18)!=(isec%18)){
2481     dalpha = -((refSector%18)-(isec%18))*TMath::TwoPi()/18.;
2482   }
2483   Double_t ca = TMath::Cos(dalpha);
2484   Double_t sa = TMath::Sin(dalpha);
2485   //
2486   //
2487   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
2488   //
2489   // full track fit parameters
2490   //
2491   static TLinearFitter fyf(2,"pol1");   // change to static - suggestion of calgrind - 30 % of time
2492   static TLinearFitter fzf(2,"pol1");   // relative to time of given class
2493   TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
2494   TMatrixD covY(4,4),covZ(4,4);
2495   Double_t chi2FacY =1;
2496   Double_t chi2FacZ =1;
2497   Int_t nf=0;
2498   //
2499   //
2500   //
2501   Float_t erry=0.1;   // initial cluster error estimate
2502   Float_t errz=0.1;   // initial cluster error estimate
2503   for (Int_t iter=0; iter<2; iter++){
2504     fyf.ClearPoints();
2505     fzf.ClearPoints();
2506     for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
2507       AliTPCclusterMI *c=track->GetClusterPointer(irow);
2508       if (!c) continue;
2509       //
2510       if (c->GetDetector()%36!=(isec%36)) continue;
2511       if (!both && c->GetDetector()!=isec) continue;
2512
2513       if (c->GetRow()<kdrow0) continue;
2514       //cluster position in reference frame
2515       Double_t lxR   =   ca*c->GetX()-sa*c->GetY();
2516       Double_t lyR   =  +sa*c->GetX()+ca*c->GetY();
2517       Double_t lzR   =  c->GetZ();
2518
2519       Double_t dx = lxR -xRef;   // distance to reference X
2520       Double_t x[2]={dx, dx*dx};
2521
2522       Double_t yfitR  =    pyf[0]+pyf[1]*dx;  // fit value Y in ref frame
2523       Double_t zfitR  =    pzf[0]+pzf[1]*dx;  // fit value Z in ref frame
2524       //
2525       Double_t yfit   =  -sa*lxR + ca*yfitR;  // fit value Y in local frame
2526       //
2527       if (iter==0 &&c->GetType()<0) continue;
2528       if (iter>0){
2529         if (TMath::Abs(lyR-yfitR)>kMaxDist*erry) continue;
2530         if (TMath::Abs(lzR-zfitR)>kMaxDist*errz) continue;
2531         Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2532         if (isec<36 && dedge<kedgeyIn) continue;
2533         if (isec>35 && dedge<kedgeyOut) continue;
2534         Double_t corrtrY =
2535           corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
2536                                   c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2537         Double_t corrclY =
2538           corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
2539                                   c->GetY(),c->GetY(), c->GetZ(), pyf[1], c->GetMax(),2.5);
2540         if (TMath::Abs((corrtrY+corrclY)*0.5)>kMaxCorrY) continue;
2541         if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2542       }
2543       fyf.AddPoint(x,lyR,erry);
2544       fzf.AddPoint(x,lzR,errz);
2545     }
2546     nf = fyf.GetNpoints();
2547     if (nf<kMinClusterF) return 0;   // not enough points - skip
2548     fyf.Eval();
2549     fyf.GetParameters(pyf);
2550     fyf.GetErrors(peyf);
2551     fzf.Eval();
2552     fzf.GetParameters(pzf);
2553     fzf.GetErrors(pezf);
2554     chi2FacY = TMath::Sqrt(fyf.GetChisquare()/(fyf.GetNpoints()-2.));
2555     chi2FacZ = TMath::Sqrt(fzf.GetChisquare()/(fzf.GetNpoints()-2.));
2556     peyf[0]*=chi2FacY;
2557     peyf[1]*=chi2FacY;
2558     pezf[0]*=chi2FacZ;
2559     pezf[1]*=chi2FacZ;
2560     erry*=chi2FacY;
2561     errz*=chi2FacZ;
2562     fyf.GetCovarianceMatrix(covY);
2563     fzf.GetCovarianceMatrix(covZ);
2564     for (Int_t i0=0;i0<2;i0++)
2565       for (Int_t i1=0;i1<2;i1++){
2566         covY(i0,i1)*=chi2FacY*chi2FacY;
2567         covZ(i0,i1)*=chi2FacZ*chi2FacZ;
2568       }
2569   }
2570   fitParam[0] = xRef;
2571   //
2572   fitParam[1] = pyf[0];
2573   fitParam[2] = pyf[1];
2574   fitParam[3] = pzf[0];
2575   fitParam[4] = pzf[1];
2576   //
2577   fitParam[5] = 0;
2578   fitParam[6] = peyf[0];
2579   fitParam[7] = peyf[1];
2580   fitParam[8] = pezf[0];
2581   fitParam[9] = pezf[1];
2582   //
2583   //
2584   tparam(0,0) = pyf[0];
2585   tparam(1,0) = pyf[1];
2586   tparam(2,0) = pzf[0];
2587   tparam(3,0) = pzf[1];
2588   //
2589   tcovar(0,0) = covY(0,0);
2590   tcovar(1,1) = covY(1,1);
2591   tcovar(1,0) = covY(1,0);
2592   tcovar(0,1) = covY(0,1);
2593   tcovar(2,2) = covZ(0,0);
2594   tcovar(3,3) = covZ(1,1);
2595   tcovar(3,2) = covZ(1,0);
2596   tcovar(2,3) = covZ(0,1);
2597   return nf;
2598 }
2599
2600 void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){
2601   //
2602   // Update the cluster residula histograms for setup with field
2603   // Kalman track fitting is used
2604   //  Only high momenta primary tracks used
2605   //
2606   // 1. Apply selection
2607   // 2. Refit the track - in-out
2608   // 3. Refit the track - out-in
2609   // 4. Combine In and Out track - - fil cluster residuals
2610   //
2611   if (!fCurrentFriendTrack) return;
2612   if (!fCurrentFriendTrack->GetTPCOut()) return;
2613   const Double_t kPtCut=1.0;    // pt
2614   const Double_t kSnpCut=0.2; // snp cut
2615   const Double_t kNclCut=120; //
2616   const Double_t kVertexCut=1;
2617   const Double_t kMaxDist=0.5; // max distance between tracks and cluster
2618   const Double_t kEdgeCut = 2.5;
2619   const Double_t kDelta2=0.2*0.2;  // initial increase in covar matrix
2620   const Double_t kSigma=0.3;       // error increase towards edges of TPC
2621   const Double_t kSkipBoundary=7.5;  // skip track updates in the boundary IFC,OFC, IO
2622   //
2623   if (!fCurrentTrack) return;
2624   if (!fCurrentFriendTrack) return;
2625   Float_t vertexXY=0,vertexZ=0;
2626   fCurrentTrack->GetImpactParameters(vertexXY,vertexZ);
2627   if (TMath::Abs(vertexXY)>kVertexCut) return;
2628   if (TMath::Abs(vertexZ)>kVertexCut) return;
2629   if (TMath::Abs(seed->Pt())<kPtCut) return;
2630   if (seed->GetNumberOfClusters()<kNclCut) return;
2631   if (TMath::Abs(seed->GetSnp())>kSnpCut) return;
2632   if (!fClusterDelta[0])  MakeResidualHistos();
2633   //
2634   AliExternalTrackParam fitIn[160];
2635   AliExternalTrackParam fitOut[160];
2636   AliTPCROC * roc = AliTPCROC::Instance();
2637   Double_t xmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
2638   Double_t xDiff     = ( -roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
2639   Double_t xIFC      = ( roc->GetPadRowRadii(0,0));
2640   Double_t xOFC      = ( roc->GetPadRowRadii(36,roc->GetNRows(36)-1));
2641   //
2642   Int_t detector=-1;
2643   //
2644   //
2645   AliExternalTrackParam trackIn  = *(fCurrentTrack->GetInnerParam());
2646   AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut());
2647   trackIn.ResetCovariance(10);
2648   trackOut.ResetCovariance(10);
2649   Double_t *covarIn = (Double_t*)trackIn.GetCovariance();
2650   Double_t *covarOut = (Double_t*)trackOut.GetCovariance();
2651   covarIn[0]+=kDelta2; covarIn[2]+=kDelta2;
2652   covarIn[5]+=kDelta2/(100.*100.); covarIn[9]=kDelta2/(100.*100.);
2653   covarIn[14]+=kDelta2/(5.*5.);
2654   covarOut[0]+=kDelta2; covarOut[2]+=kDelta2;
2655   covarOut[5]+=kDelta2/(100.*100.); covarOut[9]=kDelta2/(100.*100.);
2656   covarOut[14]+=kDelta2/(5.*5.);
2657   //
2658   static Double_t mass =    TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2659   //
2660   Int_t ncl=0;
2661   for (Int_t irow=0; irow<160; irow++){
2662     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2663     if (!cl) continue;
2664     if (cl->GetX()<80) continue;
2665     if (detector<0) detector=cl->GetDetector()%36;
2666     if (detector!=cl->GetDetector()%36) return; // cluster from different sectors
2667     // skip such tracks
2668     ncl++;
2669   }
2670   if (ncl<kNclCut) return;
2671   Int_t nclIn=0,nclOut=0;
2672   Double_t xyz[3];
2673   //
2674   // Refit out - store residual maps
2675   //
2676   for (Int_t irow=0; irow<160; irow++){
2677     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2678     if (!cl) continue;
2679     if (cl->GetX()<80) continue;
2680     if (detector<0) detector=cl->GetDetector()%36;
2681     Int_t sector = cl->GetDetector();
2682     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
2683     if (cl->GetDetector()%36!=detector) continue;
2684     if (TMath::Abs(dalpha)>0.01){
2685       if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
2686     }
2687     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
2688     Double_t cov[3]={0.1,0.,0.1};
2689     Double_t dedge =  cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackOut.GetY());
2690     Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
2691     dmiddle*=dmiddle;
2692     //
2693     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2694     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2695     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2696     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2697     cov[0]+=kSigma/dedge;      // bigger error close to the boundary
2698     cov[2]+=kSigma/dedge;      // bigger error close to the boundary
2699     cov[0]*=cov[0];
2700     cov[2]*=cov[2];
2701     if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
2702     if (TMath::Abs(dedge)<kEdgeCut) continue;
2703     //
2704     Bool_t doUpdate=kTRUE;
2705     if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
2706     if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
2707     if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
2708     //
2709     if (TMath::Abs(cl->GetY()-trackOut.GetY())<kMaxDist){
2710       nclOut++;
2711       if (doUpdate) trackOut.Update(&r[1],cov);
2712     }
2713     fitOut[irow]=trackOut;
2714   }
2715
2716   //
2717   // Refit In - store residual maps
2718   //
2719   for (Int_t irow=159; irow>=0; irow--){
2720     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2721     if (!cl) continue;
2722     if (cl->GetX()<80) continue;
2723     if (detector<0) detector=cl->GetDetector()%36;
2724     Int_t sector = cl->GetDetector();
2725     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
2726     if (cl->GetDetector()%36!=detector) continue;
2727     if (TMath::Abs(dalpha)>0.01){
2728       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
2729     }
2730     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
2731     Double_t cov[3]={0.1,0.,0.1};
2732     Double_t dedge =  cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackIn.GetY());
2733     Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
2734     dmiddle*=dmiddle;
2735     //
2736     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2737     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2738     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2739     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2740     cov[0]+=kSigma/dedge;      // bigger error close to the boundary
2741     cov[2]+=kSigma/dedge;      // bigger error close to the boundary
2742     cov[0]*=cov[0];
2743     cov[2]*=cov[2];
2744     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
2745     if (TMath::Abs(dedge)<kEdgeCut) continue;
2746     Bool_t doUpdate=kTRUE;
2747     if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
2748     if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
2749     if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
2750     if (TMath::Abs(cl->GetY()-trackIn.GetY())<kMaxDist){
2751       nclIn++;
2752       if (doUpdate) trackIn.Update(&r[1],cov);
2753     }
2754     fitIn[irow]=trackIn;
2755   }
2756   //
2757   //
2758   for (Int_t irow=159; irow>=0; irow--){
2759     //
2760     // Update kalman - +- direction
2761     // Store cluster residuals
2762     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2763     if (!cl) continue;
2764     if (cl->GetX()<80) continue;
2765     if (detector<0) detector=cl->GetDetector()%36;
2766     if (cl->GetDetector()%36!=detector) continue;
2767     if (fitIn[irow].GetX()<80) continue;
2768     if (fitOut[irow].GetX()<80) continue;
2769     AliExternalTrackParam trackSmooth = fitIn[irow];
2770     AliTrackerBase::UpdateTrack(trackSmooth, fitOut[irow]);
2771     //
2772     Double_t resVector[5];
2773     trackSmooth.GetXYZ(xyz);
2774     resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
2775     if (resVector[1]<0) resVector[1]+=18;
2776     resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
2777     resVector[3]= cl->GetZ()/resVector[2];
2778     //
2779     resVector[0]= cl->GetY()-trackSmooth.GetY();
2780     fClusterDelta[0]->Fill(resVector);
2781     resVector[0]= cl->GetZ()-trackSmooth.GetZ();
2782     fClusterDelta[1]->Fill(resVector);
2783   }
2784
2785 }
2786
2787
2788 void  AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
2789   //
2790   // Update Kalman filter of Alignment - only setup without filed
2791   //       IROC - OROC quadrants
2792   //
2793   if (TMath::Abs(AliTracker::GetBz())>0.5) return;
2794   if (!fClusterDelta[0])  MakeResidualHistos();
2795   //  const Int_t kMinClusterF=40;
2796   const Int_t kMinClusterFit=10;
2797   const Int_t kMinClusterQ=10;
2798   //
2799   const  Int_t     kdrow1Fit =5;         // rows to skip from fit at the end
2800   const  Int_t     kdrow0Fit =10;        // rows to skip from fit at beginning
2801   const  Float_t   kedgey=3.0;
2802   const  Float_t   kMaxDist=1;
2803   const  Float_t   kMaxCorrY=0.05;
2804   const  Float_t   kPRFWidth = 0.6;   //cut 2 sigma of PRF
2805   isec = isec%36;     // use the hardware numbering
2806   //
2807   //
2808   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
2809   //
2810   // full track fit parameters
2811   //
2812   static TLinearFitter fyf(2,"pol1");   // make it static - too much time for comiling of formula
2813   static TLinearFitter fzf(2,"pol1");   // calgrind recomendation
2814   TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
2815   TVectorD pyfc(2),pzfc(2);
2816   Int_t nf=0;
2817   //
2818   // make full fit as reference
2819   //
2820   for (Int_t iter=0; iter<2; iter++){
2821     fyf.ClearPoints();
2822     fzf.ClearPoints();
2823     for (Int_t irow=kdrow0Fit;irow<159-kdrow1Fit;irow++) {
2824       AliTPCclusterMI *c=track->GetClusterPointer(irow);
2825       if (!c) continue;
2826       if ((c->GetDetector()%36)!=isec) continue;
2827       if (c->GetRow()<kdrow0Fit) continue;
2828       Double_t dx = c->GetX()-fXmiddle;
2829       Double_t x[2]={dx, dx*dx};
2830       if (iter==0 &&c->GetType()<0) continue;
2831       if (iter==1){
2832         Double_t yfit  =  pyf[0]+pyf[1]*dx;
2833         Double_t zfit  =  pzf[0]+pzf[1]*dx;
2834         Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2835         if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
2836         if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
2837         if (dedge<kedgey) continue;
2838         Double_t corrtrY =
2839           corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
2840                                   c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2841         if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2842       }
2843       if (TMath::Abs(x[0])<10){
2844         fyf.AddPoint(x,c->GetY(),0.1); //use only middle rows+-10cm
2845         fzf.AddPoint(x,c->GetZ(),0.1);
2846       }
2847     }
2848     nf = fyf.GetNpoints();
2849     if (fyf.GetNpoints()<kMinClusterFit) return;   // not enough points - skip
2850     if (fzf.GetNpoints()<kMinClusterFit) return;   // not enough points - skip
2851     fyf.Eval();
2852     fyf.GetParameters(pyf);
2853     fyf.GetErrors(peyf);
2854     fzf.Eval();
2855     fzf.GetParameters(pzf);
2856     fzf.GetErrors(pezf);
2857   }
2858   //
2859   //
2860   //
2861   TVectorD vecX(160);          // x         vector
2862   TVectorD vecY(160);          // residuals vector
2863   TVectorD vecZ(160);                              // residuals vector
2864   TVectorD vPosG(3);                  //vertex position
2865   TVectorD vPosL(3);                 // vertex position in the TPC local system
2866   TVectorF vImpact(2);               //track impact parameter
2867   //  Double_t tofSignal= fCurrentTrack->GetTOFsignal();      // tof signal
2868   TVectorF tpcPosG(3);                                    // global position of track at the middle of fXmiddle
2869   Double_t lphi = TMath::ATan2(pyf[0],fXmiddle);          // expected local phi angle - if vertex at 0
2870   Double_t gphi = 2.*TMath::Pi()*(isec%18+0.5)/18.+lphi;  // expected global phi if vertex at 0
2871   Double_t ky   = pyf[0]/fXmiddle;
2872   Double_t kyE  =0, kzE=0;    // ky and kz expected
2873   Double_t alpha =2.*TMath::Pi()*(isec%18+0.5)/18.;
2874   Double_t scos=TMath::Cos(alpha);
2875   Double_t ssin=TMath::Sin(alpha);
2876   const AliESDVertex* vertex = fCurrentEvent->GetPrimaryVertexTracks();
2877   vertex->GetXYZ(vPosG.GetMatrixArray());
2878   fCurrentTrack->GetImpactParameters(vImpact[0],vImpact[1]);  // track impact parameters
2879   //
2880   tpcPosG[0]= scos*fXmiddle-ssin*pyf[0];
2881   tpcPosG[1]=+ssin*fXmiddle+scos*pyf[0];
2882   tpcPosG[2]=pzf[0];
2883   vPosL[0]= scos*vPosG[0]+ssin*vPosG[1];
2884   vPosL[1]=-ssin*vPosG[0]+scos*vPosG[1];
2885   vPosL[2]= vPosG[2];
2886   kyE = (pyf[0]-vPosL[1])/(fXmiddle-vPosL[0]);
2887   kzE = (pzf[0]-vPosL[2])/(fXmiddle-vPosL[0]);
2888   //
2889   // get constrained parameters
2890   //
2891   Double_t xvertex=vPosL[0]-fXmiddle;
2892   fyf.AddPoint(&xvertex,vPosL[1], 0.00001);
2893   fzf.AddPoint(&xvertex,vPosL[2], 2.);
2894   fyf.Eval();
2895   fyf.GetParameters(pyfc);
2896   fzf.Eval();
2897   fzf.GetParameters(pzfc);
2898   //
2899   //
2900   // Make Fitters and params for 5 fitters
2901   // 1-4 OROC quadrants
2902   //   0 IROC
2903   //
2904   static TLinearFitter *fittersY[5]={0,0,0,0,0};   // calgrind recomendation - fater to clear points
2905   static TLinearFitter *fittersZ[5]={0,0,0,0,0};   // than create the fitter
2906   if (fittersY[0]==0){
2907     for (Int_t i=0;i<5;i++) {
2908       fittersY[i] = new TLinearFitter(2,"pol1");
2909       fittersZ[i] = new TLinearFitter(2,"pol1");
2910     }
2911   }
2912   //
2913   Int_t npoints[5];
2914   TVectorD paramsY[5];
2915   TVectorD errorsY[5];
2916   TMatrixD covY[5];
2917   Double_t chi2CY[5];
2918   TVectorD paramsZ[5];
2919   TVectorD errorsZ[5];
2920   TMatrixD covZ[5];
2921   Double_t chi2CZ[5];
2922   for (Int_t i=0;i<5;i++) {
2923     npoints[i]=0;
2924     paramsY[i].ResizeTo(2);
2925     errorsY[i].ResizeTo(2);
2926     covY[i].ResizeTo(2,2);
2927     paramsZ[i].ResizeTo(2);
2928     errorsZ[i].ResizeTo(2);
2929     covZ[i].ResizeTo(2,2);
2930     fittersY[i]->ClearPoints();
2931     fittersZ[i]->ClearPoints();
2932   }
2933   //
2934   // Update fitters
2935   //
2936   Int_t countRes=0;
2937   for (Int_t irow=0;irow<159;irow++) {
2938     AliTPCclusterMI *c=track->GetClusterPointer(irow);
2939     if (!c) continue;
2940     if ((c->GetDetector()%36)!=isec) continue;
2941     Double_t dx = c->GetX()-fXmiddle;
2942     Double_t x[2]={dx, dx*dx};
2943     Double_t yfit  =  pyf[0]+pyf[1]*dx;
2944     Double_t zfit  =  pzf[0]+pzf[1]*dx;
2945     Double_t yfitC  =  pyfc[0]+pyfc[1]*dx;
2946     Double_t zfitC  =  pzfc[0]+pzfc[1]*dx;
2947     Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2948     if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
2949     if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
2950     if (dedge<kedgey) continue;
2951     Double_t corrtrY =
2952       corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
2953                               c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2954     if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2955     //
2956     vecX[countRes]=c->GetX();
2957     vecY[countRes]=c->GetY()-yfit;
2958     vecZ[countRes]=c->GetZ()-zfit;
2959     countRes++;
2960     //
2961     // Fill THnSparse cluster residuals
2962     // use only primary candidates with ITS signal
2963     if (fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){
2964       Double_t resVector[5];
2965       resVector[1]= 9.*gphi/TMath::Pi();
2966       resVector[2]= TMath::Sqrt(c->GetX()*c->GetX()+c->GetY()*c->GetY());
2967       resVector[3]= c->GetZ()/resVector[2];
2968       //
2969       //
2970       resVector[0]= c->GetY()-yfitC;
2971       fClusterDelta[0]->Fill(resVector);
2972       resVector[0]= c->GetZ()-zfitC;
2973       fClusterDelta[1]->Fill(resVector);
2974     }
2975     if (c->GetRow()<kdrow0Fit) continue;
2976     if (c->GetRow()>159-kdrow1Fit) continue;
2977     //
2978
2979     if (c->GetDetector()>35){
2980       if (c->GetX()<fXquadrant){
2981         if (yfit<-kPRFWidth) fittersY[1]->AddPoint(x,c->GetY(),0.1);
2982         if (yfit<-kPRFWidth) fittersZ[1]->AddPoint(x,c->GetZ(),0.1);
2983         if (yfit>kPRFWidth) fittersY[2]->AddPoint(x,c->GetY(),0.1);
2984         if (yfit>kPRFWidth) fittersZ[2]->AddPoint(x,c->GetZ(),0.1);
2985       }
2986       if (c->GetX()>fXquadrant){
2987         if (yfit<-kPRFWidth) fittersY[3]->AddPoint(x,c->GetY(),0.1);
2988         if (yfit<-kPRFWidth) fittersZ[3]->AddPoint(x,c->GetZ(),0.1);
2989         if (yfit>kPRFWidth) fittersY[4]->AddPoint(x,c->GetY(),0.1);
2990         if (yfit>kPRFWidth) fittersZ[4]->AddPoint(x,c->GetZ(),0.1);
2991       }
2992     }
2993     if (c->GetDetector()<36){
2994       fittersY[0]->AddPoint(x,c->GetY(),0.1);
2995       fittersZ[0]->AddPoint(x,c->GetZ(),0.1);
2996     }
2997   }
2998   //
2999   //
3000   //
3001   for (Int_t i=0;i<5;i++) {
3002     npoints[i] = fittersY[i]->GetNpoints();
3003     if (npoints[i]>=kMinClusterQ){
3004       // Y fit
3005       fittersY[i]->Eval();
3006       Double_t chi2FacY = TMath::Sqrt(fittersY[i]->GetChisquare()/(fittersY[i]->GetNpoints()-2));
3007       chi2CY[i]=chi2FacY;
3008       fittersY[i]->GetParameters(paramsY[i]);
3009       fittersY[i]->GetErrors(errorsY[i]);
3010       fittersY[i]->GetCovarianceMatrix(covY[i]);
3011       // renormalize errors
3012       errorsY[i][0]*=chi2FacY;
3013       errorsY[i][1]*=chi2FacY;
3014       covY[i](0,0)*=chi2FacY*chi2FacY;
3015       covY[i](0,1)*=chi2FacY*chi2FacY;
3016       covY[i](1,0)*=chi2FacY*chi2FacY;
3017       covY[i](1,1)*=chi2FacY*chi2FacY;
3018       // Z fit
3019       fittersZ[i]->Eval();
3020       Double_t chi2FacZ = TMath::Sqrt(fittersZ[i]->GetChisquare()/(fittersZ[i]->GetNpoints()-2));
3021       chi2CZ[i]=chi2FacZ;
3022       fittersZ[i]->GetParameters(paramsZ[i]);
3023       fittersZ[i]->GetErrors(errorsZ[i]);
3024       fittersZ[i]->GetCovarianceMatrix(covZ[i]);
3025       // renormalize errors
3026       errorsZ[i][0]*=chi2FacZ;
3027       errorsZ[i][1]*=chi2FacZ;
3028       covZ[i](0,0)*=chi2FacZ*chi2FacZ;
3029       covZ[i](0,1)*=chi2FacZ*chi2FacZ;
3030       covZ[i](1,0)*=chi2FacZ*chi2FacZ;
3031       covZ[i](1,1)*=chi2FacZ*chi2FacZ;
3032     }
3033   }
3034   //
3035   // void UpdateSectorKalman
3036   //
3037   for (Int_t i0=0;i0<5;i0++){
3038     for (Int_t i1=i0+1;i1<5;i1++){
3039       if(npoints[i0]<kMinClusterQ) continue;
3040       if(npoints[i1]<kMinClusterQ) continue;
3041       TMatrixD v0(4,1),v1(4,1);        // measurement
3042       TMatrixD cov0(4,4),cov1(4,4);    // covariance
3043       //
3044       v0(0,0)= paramsY[i0][0];         v1(0,0)= paramsY[i1][0];
3045       v0(1,0)= paramsY[i0][1];         v1(1,0)= paramsY[i1][1];
3046       v0(2,0)= paramsZ[i0][0];         v1(2,0)= paramsZ[i1][0];
3047       v0(3,0)= paramsZ[i0][1];         v1(3,0)= paramsZ[i1][1];
3048       //covariance i0
3049       cov0(0,0) = covY[i0](0,0);
3050       cov0(1,1) = covY[i0](1,1);
3051       cov0(1,0) = covY[i0](1,0);
3052       cov0(0,1) = covY[i0](0,1);
3053       cov0(2,2) = covZ[i0](0,0);
3054       cov0(3,3) = covZ[i0](1,1);
3055       cov0(3,2) = covZ[i0](1,0);
3056       cov0(2,3) = covZ[i0](0,1);
3057       //covariance i1
3058       cov1(0,0) = covY[i1](0,0);
3059       cov1(1,1) = covY[i1](1,1);
3060       cov1(1,0) = covY[i1](1,0);
3061       cov1(0,1) = covY[i1](0,1);
3062       cov1(2,2) = covZ[i1](0,0);
3063       cov1(3,3) = covZ[i1](1,1);
3064       cov1(3,2) = covZ[i1](1,0);
3065       cov1(2,3) = covZ[i1](0,1);
3066       //
3067       // And now update
3068       //
3069       if (TMath::Abs(pyf[1])<0.8){ //angular cut
3070         UpdateSectorKalman(isec, i0,i1, &v0,&cov0,&v1,&cov1);
3071       }
3072     }
3073   }
3074
3075   //
3076   // Dump debug information
3077   //
3078   if (fStreamLevel>0){
3079     // get vertex position
3080      //
3081     TTreeSRedirector *cstream = GetDebugStreamer();
3082
3083
3084     if (cstream){
3085       for (Int_t i0=0;i0<5;i0++){
3086         for (Int_t i1=i0;i1<5;i1++){
3087           if (i0==i1) continue;
3088           if(npoints[i0]<kMinClusterQ) continue;
3089           if(npoints[i1]<kMinClusterQ) continue;
3090           (*cstream)<<"sectorAlign"<<
3091             "run="<<fRun<<              //  run number
3092             "event="<<fEvent<<          //  event number
3093             "time="<<fTime<<            //  time stamp of event
3094             "trigger="<<fTrigger<<      //  trigger
3095             "triggerClass="<<&fTriggerClass<<      //  trigger
3096             "mag="<<fMagF<<             //  magnetic field
3097             "isec="<<isec<<             //  current sector
3098             //
3099             "xref="<<fXmiddle<<         // reference X
3100             "vPosG.="<<&vPosG<<        // vertex position in global system
3101             "vPosL.="<<&vPosL<<        // vertex position in local  system
3102             "vImpact.="<<&vImpact<<   // track impact parameter
3103             //"tofSignal="<<tofSignal<<   // tof signal
3104             "tpcPosG.="<<&tpcPosG<<   // global position of track at the middle of fXmiddle
3105             "lphi="<<lphi<<             // expected local phi angle - if vertex at 0
3106             "gphi="<<gphi<<             // expected global phi if vertex at 0
3107             "ky="<<ky<<
3108             "kyE="<<kyE<<               // expect ky - assiming pirmary track
3109             "kzE="<<kzE<<               // expected kz - assuming primary tracks
3110             "salpha="<<alpha<<          // sector alpha
3111             "scos="<<scos<<              // tracking cosinus
3112             "ssin="<<ssin<<              // tracking sinus
3113             //
3114             // full fit
3115             //
3116             "nf="<<nf<<                 //  total number of points
3117             "pyf.="<<&pyf<<             //  full OROC fit y
3118             "pzf.="<<&pzf<<             //  full OROC fit z
3119             "vX.="<<&vecX<<              //  x cluster
3120             "vY.="<<&vecY<<              //  y residual cluster
3121             "vZ.="<<&vecZ<<              //  z residual cluster
3122             // quadrant and IROC fit
3123             "i0="<<i0<<                 // quadrant number
3124             "i1="<<i1<<
3125             "n0="<<npoints[i0]<<        // number of points
3126             "n1="<<npoints[i1]<<
3127             //
3128             "py0.="<<&paramsY[i0]<<       // parameters
3129             "py1.="<<&paramsY[i1]<<
3130             "ey0.="<<&errorsY[i0]<<       // errors
3131             "ey1.="<<&errorsY[i1]<<
3132             "chiy0="<<chi2CY[i0]<<       // chi2s
3133             "chiy1="<<chi2CY[i1]<<
3134             //
3135             "pz0.="<<&paramsZ[i0]<<       // parameters
3136             "pz1.="<<&paramsZ[i1]<<
3137             "ez0.="<<&errorsZ[i0]<<       // errors
3138             "ez1.="<<&errorsZ[i1]<<
3139             "chiz0="<<chi2CZ[i0]<<       // chi2s
3140             "chiz1="<<chi2CZ[i1]<<
3141             "\n";
3142         }
3143       }
3144     }
3145   }
3146 }
3147
3148 void AliTPCcalibAlign::UpdateSectorKalman(Int_t sector, Int_t quadrant0, Int_t quadrant1,  TMatrixD *const p0, TMatrixD *const c0, TMatrixD *const p1, TMatrixD *const c1 ){
3149   //
3150   //
3151   // tracks are refitted at sector middle
3152   //
3153   if (fArraySectorIntParam.At(0)==NULL) MakeSectorKalman();
3154   //
3155   //
3156   static TMatrixD matHk(4,30);    // vector to mesurement
3157   static TMatrixD measR(4,4);     // measurement error
3158   //  static TMatrixD matQk(2,2);     // prediction noise vector
3159   static TMatrixD vecZk(4,1);     // measurement
3160   //
3161   static TMatrixD vecYk(4,1);     // Innovation or measurement residual
3162   static TMatrixD matHkT(30,4);   // helper matrix Hk transpose
3163   static TMatrixD matSk(4,4);     // Innovation (or residual) covariance
3164   static TMatrixD matKk(30,4);    // Optimal Kalman gain
3165   static TMatrixD mat1(30,30);    // update covariance matrix
3166   static TMatrixD covXk2(30,30);  // helper matrix
3167   //
3168   TMatrixD *vOrig = (TMatrixD*)(fArraySectorIntParam.At(sector));
3169   TMatrixD *cOrig = (TMatrixD*)(fArraySectorIntCovar.At(sector));
3170   //
3171   TMatrixD vecXk(*vOrig);             // X vector
3172   TMatrixD covXk(*cOrig);             // X covariance
3173   //
3174   //Unit matrix
3175   //
3176   for (Int_t i=0;i<30;i++)
3177     for (Int_t j=0;j<30;j++){
3178       mat1(i,j)=0;
3179       if (i==j) mat1(i,j)=1;
3180     }
3181   //
3182   //
3183   // matHk - vector to measurement
3184   //
3185   for (Int_t i=0;i<4;i++)
3186     for (Int_t j=0;j<30;j++){
3187       matHk(i,j)=0;
3188     }
3189   //
3190   // Measurement
3191   // 0  - y
3192   // 1  - ky
3193   // 2  - z
3194   // 3  - kz
3195
3196   matHk(0,6*quadrant1+4)  =  1.;            // delta y
3197   matHk(1,6*quadrant1+0)  =  1.;            // delta ky
3198   matHk(2,6*quadrant1+5)  =  1.;            // delta z
3199   matHk(3,6*quadrant1+1)  =  1.;            // delta kz
3200   // bug fix 24.02  - aware of sign in dx
3201   matHk(0,6*quadrant1+3)  =  -(*p0)(1,0);    // delta x to delta y  - through ky
3202   matHk(2,6*quadrant1+3)  =  -(*p0)(3,0);    // delta x to delta z  - thorugh kz
3203   matHk(2,6*quadrant1+2)  =  ((*p0)(0,0));  // y       to delta z  - through psiz
3204   //
3205   matHk(0,6*quadrant0+4)  =  -1.;           // delta y
3206   matHk(1,6*quadrant0+0)  =  -1.;           // delta ky
3207   matHk(2,6*quadrant0+5)  =  -1.;           // delta z
3208   matHk(3,6*quadrant0+1)  =  -1.;           // delta kz
3209   // bug fix 24.02 be aware of sign in dx
3210   matHk(0,6*quadrant0+3)  =  ((*p0)(1,0)); // delta x to delta y  - through ky
3211   matHk(2,6*quadrant0+3)  =  ((*p0)(3,0)); // delta x to delta z  - thorugh kz
3212   matHk(2,6*quadrant0+2)  =  -((*p0)(0,0)); // y       to delta z  - through psiz
3213
3214   //
3215   //
3216
3217   vecZk =(*p1)-(*p0);                 // measurement
3218   measR =(*c1)+(*c0);                 // error of measurement
3219   vecYk = vecZk-matHk*vecXk;          // Innovation or measurement residual
3220   //
3221   //
3222   matHkT=matHk.T(); matHk.T();
3223   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
3224   matSk.Invert();
3225   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
3226   vecXk += matKk*vecYk;                    //  updated vector
3227   covXk2= (mat1-(matKk*matHk));
3228   covXk =  covXk2*covXk;
3229   //
3230   //
3231   (*cOrig)=covXk;
3232   (*vOrig)=vecXk;
3233 }
3234
3235 void AliTPCcalibAlign::MakeSectorKalman(){
3236   //
3237   // Make a initial Kalman paramaters for IROC - Quadrants alignment
3238   //
3239   TMatrixD param(5*6,1);
3240   TMatrixD covar(5*6,5*6);
3241   //
3242   // Set inital parameters
3243   //
3244   for (Int_t ip=0;ip<5*6;ip++) param(ip,0)=0;  // mean alignment to 0
3245   //
3246   for (Int_t iq=0;iq<5;iq++){
3247     // Initial uncertinty
3248     covar(iq*6+0,iq*6+0) = 0.002*0.002;        // 2 mrad
3249     covar(iq*6+1,iq*6+1) = 0.002*0.002;        // 2 mrad  rotation
3250     covar(iq*6+2,iq*6+2) = 0.002*0.002;        // 2 mrad
3251     //
3252     covar(iq*6+3,iq*6+3) = 0.02*0.02;        // 0.2 mm
3253     covar(iq*6+4,iq*6+4) = 0.02*0.02;        // 0.2 mm  translation
3254     covar(iq*6+5,iq*6+5) = 0.02*0.02;        // 0.2 mm
3255   }
3256
3257   for (Int_t isec=0;isec<36;isec++){
3258     fArraySectorIntParam.AddAt(param.Clone(),isec);
3259     fArraySectorIntCovar.AddAt(covar.Clone(),isec);
3260   }
3261 }
3262
3263 void AliTPCcalibAlign::UpdateSectorKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){
3264   //
3265   // Update kalman vector para0 with vector par1
3266   // Used for merging
3267   //
3268   static TMatrixD matHk(30,30);    // vector to mesurement
3269   static TMatrixD measR(30,30);     // measurement error
3270   static TMatrixD vecZk(30,1);     // measurement
3271   //
3272   static TMatrixD vecYk(30,1);     // Innovation or measurement residual
3273   static TMatrixD matHkT(30,30);   // helper matrix Hk transpose
3274   static TMatrixD matSk(30,30);     // Innovation (or residual) covariance
3275   static TMatrixD matKk(30,30);    // Optimal Kalman gain
3276   static TMatrixD mat1(30,30);    // update covariance matrix
3277   static TMatrixD covXk2(30,30);  // helper matrix
3278   //
3279   TMatrixD vecXk(par0);             // X vector
3280   TMatrixD covXk(cov0);             // X covariance
3281
3282   //
3283   //Unit matrix
3284   //
3285   for (Int_t i=0;i<30;i++)
3286     for (Int_t j=0;j<30;j++){
3287       mat1(i,j)=0;
3288       if (i==j) mat1(i,j)=1;
3289     }
3290   matHk = mat1;                       // unit matrix
3291   //
3292   vecZk = par1;                       // measurement
3293   measR = cov1;                        // error of measurement
3294   vecYk = vecZk-matHk*vecXk;          // Innovation or measurement residual
3295   //
3296   matHkT=matHk.T(); matHk.T();
3297   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
3298   matSk.Invert();
3299   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
3300   //matKk.Print();
3301   vecXk += matKk*vecYk;                    //  updated vector
3302   covXk2= (mat1-(matKk*matHk));
3303   covXk =  covXk2*covXk;
3304   CheckCovariance(covXk);
3305   CheckCovariance(cov1);
3306   //
3307   par0  = vecXk;                      // update measurement param
3308   cov0  = covXk;                      // update measurement covar
3309 }
3310
3311 Double_t AliTPCcalibAlign::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
3312   //
3313   // Get position correction for given sector
3314   //
3315
3316   TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
3317   if (!param) return 0;
3318   Int_t quadrant=0;
3319   if(lx>fXIO){
3320     if (lx<fXquadrant) {
3321       if (ly<0) quadrant=1;
3322       if (ly>0) quadrant=2;
3323     }
3324     if (lx>fXquadrant) {
3325       if (ly<0) quadrant=3;
3326       if (ly>0) quadrant=4;
3327     }
3328   }
3329   Double_t a10 = (*param)(quadrant*6+0,0);
3330   Double_t a20 = (*param)(quadrant*6+1,0);
3331   Double_t a21 = (*param)(quadrant*6+2,0);
3332   Double_t dx  = (*param)(quadrant*6+3,0);
3333   Double_t dy  = (*param)(quadrant*6+4,0);
3334   Double_t dz  = (*param)(quadrant*6+5,0);
3335   Double_t deltaX = lx-fXIO;
3336   if (coord==0) return dx;
3337   if (coord==1) return dy+deltaX*a10;
3338   if (coord==2) return dz+deltaX*a20+ly*a21;
3339   return 0;
3340 }
3341
3342 Double_t AliTPCcalibAlign::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
3343   //
3344   //
3345   //
3346   if (!Instance()) return 0;
3347   return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz);
3348 }
3349
3350 void AliTPCcalibAlign::MakeKalman(){
3351   //
3352   // Make a initial Kalman paramaters for sector Alignemnt
3353   //
3354   fSectorParamA = new TMatrixD(6*36+6,1);
3355   fSectorParamC = new TMatrixD(6*36+6,1);
3356   fSectorCovarA = new TMatrixD(6*36+6,6*36+6);
3357   fSectorCovarC = new TMatrixD(6*36+6,6*36+6);
3358   //
3359   // set starting parameters at 0
3360   //
3361   for (Int_t isec=0;isec<37;isec++)
3362     for (Int_t ipar=0;ipar<6;ipar++){
3363       (*fSectorParamA)(isec*6+ipar,0) =0;
3364       (*fSectorParamC)(isec*6+ipar,0) =0;
3365   }
3366   //
3367   // set starting covariance
3368   //
3369   for (Int_t isec=0;isec<36;isec++)
3370     for (Int_t ipar=0;ipar<6;ipar++){
3371       if (ipar<3){
3372         (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.002*0.002;   // 2 mrad
3373         (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.002*0.002;
3374       }
3375       if (ipar>=3){
3376         (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.02*0.02;     // 0.2 mm
3377         (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.02*0.02;
3378       }
3379   }
3380   (*fSectorCovarA)(36*6+0,36*6+0) =0.04;   // common shift y  up-up
3381   (*fSectorCovarA)(36*6+1,36*6+1) =0.04;   // common shift y  down-down
3382   (*fSectorCovarA)(36*6+2,36*6+2) =0.04;   // common shift y  up-down
3383   (*fSectorCovarA)(36*6+3,36*6+3) =0.004;   // common shift phi  up-up
3384   (*fSectorCovarA)(36*6+4,36*6+4) =0.004;   // common shift phi down-down
3385   (*fSectorCovarA)(36*6+5,36*6+5) =0.004;   // common shift phi up-down
3386   //
3387   (*fSectorCovarC)(36*6+0,36*6+0) =0.04;   // common shift y  up-up
3388   (*fSectorCovarC)(36*6+1,36*6+1) =0.04;   // common shift y  down-down
3389   (*fSectorCovarC)(36*6+2,36*6+2) =0.04;   // common shift y  up-down
3390   (*fSectorCovarC)(36*6+3,36*6+3) =0.004;   // common shift phi  up-up
3391   (*fSectorCovarC)(36*6+4,36*6+4) =0.004;   // common shift phi down-down
3392   (*fSectorCovarC)(36*6+5,36*6+5) =0.004;   // common shift phi up-down
3393 }
3394
3395 void AliTPCcalibAlign::UpdateKalman(Int_t sector0, Int_t sector1,  TMatrixD &p0, TMatrixD &c0, TMatrixD &p1, TMatrixD &c1){
3396   //
3397   // Update Kalman parameters
3398   // Note numbering from 0..36  0..17 IROC 18..35 OROC
3399   //
3400   //
3401   if (fSectorParamA==NULL) MakeKalman();
3402   if (CheckCovariance(c0)>0) return;
3403   if (CheckCovariance(c1)>0) return;
3404   const Int_t nelem = 6*36+6;
3405   //
3406   //
3407   static TMatrixD matHk(4,nelem);    // vector to mesurement
3408   static TMatrixD measR(4,4);     // measurement error
3409   static TMatrixD vecZk(4,1);     // measurement
3410   //
3411   static TMatrixD vecYk(4,1);     // Innovation or measurement residual
3412   static TMatrixD matHkT(nelem,4);   // helper matrix Hk transpose
3413   static TMatrixD matSk(4,4);     // Innovation (or residual) covariance
3414   static TMatrixD matKk(nelem,4);    // Optimal Kalman gain
3415   static TMatrixD mat1(nelem,nelem);    // update covariance matrix
3416   static TMatrixD covXk2(nelem,nelem);  // helper matrix
3417   //
3418   TMatrixD *vOrig = 0;
3419   TMatrixD *cOrig = 0;
3420   vOrig = (sector0%36>=18) ? fSectorParamA:fSectorParamC;
3421   cOrig = (sector0%36>=18) ? fSectorCovarA:fSectorCovarC;
3422   //
3423   Int_t sec0= sector0%18;
3424   Int_t sec1= sector1%18;
3425   if (sector0>35) sec0+=18;
3426   if (sector1>35) sec1+=18;
3427   //
3428   TMatrixD vecXk(*vOrig);             // X vector
3429   TMatrixD covXk(*cOrig);             // X covariance
3430   //
3431   //Unit matrix
3432   //
3433   for (Int_t i=0;i<nelem;i++)
3434     for (Int_t j=0;j<nelem;j++){
3435       mat1(i,j)=0;
3436       if (i==j) mat1(i,j)=1;
3437     }
3438   //
3439   //
3440   // matHk - vector to measurement
3441   //
3442   for (Int_t i=0;i<4;i++)
3443     for (Int_t j=0;j<nelem;j++){
3444       matHk(i,j)=0;
3445     }
3446   //
3447   // Measurement
3448   // 0  - y
3449   // 1  - ky
3450   // 2  - z
3451   // 3  - kz
3452
3453   matHk(0,6*sec1+4)  =  1.;            // delta y
3454   matHk(1,6*sec1+0)  =  1.;            // delta ky
3455   matHk(2,6*sec1+5)  =  1.;            // delta z
3456   matHk(3,6*sec1+1)  =  1.;            // delta kz
3457   matHk(0,6*sec1+3)  =  p0(1,0);    // delta x to delta y  - through ky
3458   matHk(2,6*sec1+3)  =  p0(3,0);    // delta x to delta z  - thorugh kz
3459   matHk(2,6*sec1+2)  =  p0(0,0);    // y       to delta z  - through psiz
3460   //
3461   matHk(0,6*sec0+4)  =  -1.;           // delta y
3462   matHk(1,6*sec0+0)  =  -1.;           // delta ky
3463   matHk(2,6*sec0+5)  =  -1.;           // delta z
3464   matHk(3,6*sec0+1)  =  -1.;           // delta kz
3465   matHk(0,6*sec0+3)  =  -p0(1,0); // delta x to delta y  - through ky
3466   matHk(2,6*sec0+3)  =  -p0(3,0); // delta x to delta z  - thorugh kz
3467   matHk(2,6*sec0+2)  =  -p0(0,0); // y       to delta z  - through psiz
3468
3469   Int_t dsec = (sector1%18)-(sector0%18);
3470   if (dsec<-2) dsec+=18;
3471   if (TMath::Abs(dsec)==1){
3472     //
3473     // Left right systematic fit part
3474     //
3475     Double_t dir = 0;
3476     if (dsec>0) dir= 1.;
3477     if (dsec<0) dir=-1.;
3478     if (sector0>35&&sector1>35){
3479       matHk(0,36*6+0)=dir;
3480       matHk(1,36*6+3+0)=dir;
3481     }
3482     if (sector0<36&&sector1<36){
3483       matHk(0,36*6+1)=dir;
3484       matHk(1,36*6+3+1)=dir;
3485     }
3486     if (sector0<36&&sector1>35){
3487       matHk(0,36*6+2)=dir;
3488       matHk(1,36*6+3+2)=dir;
3489     }
3490     if (sector0>35&&sector1<36){
3491       matHk(0,36*6+2)=-dir;
3492       matHk(1,36*6+3+2)=-dir;
3493     }
3494   }
3495   //
3496   //
3497   vecZk =(p1)-(p0);                 // measurement
3498   measR =(c1)+(c0);                 // error of measurement
3499   vecYk = vecZk-matHk*vecXk;          // Innovation or measurement residual
3500   //
3501   //
3502   matHkT=matHk.T(); matHk.T();
3503   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
3504   matSk.Invert();
3505   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
3506   vecXk += matKk*vecYk;                    //  updated vector
3507   covXk2= (mat1-(matKk*matHk));
3508   covXk =  covXk2*covXk;
3509
3510   if (CheckCovariance(covXk)>0) return;
3511
3512   //
3513   //
3514   (*cOrig)=covXk;
3515   (*vOrig)=vecXk;
3516 }
3517
3518
3519 void AliTPCcalibAlign::UpdateKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){
3520   //
3521   // Update kalman vector para0 with vector par1
3522   // Used for merging
3523   //
3524   Int_t nelem = 6*36+6;
3525   static TMatrixD matHk(nelem,nelem);    // vector to mesurement
3526   static TMatrixD measR(nelem,nelem);     // measurement error
3527   static TMatrixD vecZk(nelem,1);     // measurement
3528   //
3529   static TMatrixD vecYk(nelem,1);     // Innovation or measurement residual
3530   static TMatrixD matHkT(nelem,nelem);   // helper matrix Hk transpose
3531   static TMatrixD matSk(nelem,nelem);     // Innovation (or residual) covariance
3532   static TMatrixD matKk(nelem,nelem);    // Optimal Kalman gain
3533   static TMatrixD mat1(nelem,nelem);    // update covariance matrix
3534   static TMatrixD covXk2(nelem,nelem);  // helper matrix
3535   //
3536   TMatrixD vecXk(par0);             // X vector
3537   TMatrixD covXk(cov0);             // X covariance
3538
3539   //
3540   //Unit matrix
3541   //
3542   for (Int_t i=0;i<nelem;i++)
3543     for (Int_t j=0;j<nelem;j++){
3544       mat1(i,j)=0;
3545       if (i==j) mat1(i,j)=1;
3546     }
3547   matHk = mat1;                       // unit matrix
3548   //
3549   vecZk = par1;                       // measurement
3550   measR = cov1;                        // error of measurement
3551   vecYk = vecZk-matHk*vecXk;          // Innovation or measurement residual
3552   //
3553   matHkT=matHk.T(); matHk.T();
3554   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
3555   matSk.Invert();
3556   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
3557   //matKk.Print();
3558   vecXk += matKk*vecYk;                    //  updated vector
3559   covXk2= (mat1-(matKk*matHk));
3560   covXk =  covXk2*covXk;
3561   //
3562   CheckCovariance(cov0);
3563   CheckCovariance(cov1);
3564   CheckCovariance(covXk);
3565   //
3566   par0  = vecXk;                      // update measurement param
3567   cov0  = covXk;                      // update measurement covar
3568 }
3569
3570
3571
3572
3573 Int_t  AliTPCcalibAlign::CheckCovariance(TMatrixD &covar){
3574   //
3575   // check the consistency of covariance matrix
3576   //
3577   Int_t ncols = covar.GetNcols();
3578   Int_t nrows= covar.GetNrows();
3579   const Float_t kEpsilon = 0.0001;
3580   Int_t nerrors =0;
3581   //
3582   //
3583   //
3584   if (nrows!=ncols) {
3585     printf("Error 0 - wrong matrix\n");
3586     nerrors++;
3587   }
3588   //
3589   // 1. Check that the non diagonal elements
3590   //
3591   for (Int_t i0=0;i0<nrows;i0++)
3592     for (Int_t i1=i0+1;i1<ncols;i1++){
3593       Double_t r0 = covar(i0,i1)/TMath::Sqrt(covar(i0,i0)*covar(i1,i1));
3594       Double_t r1 = covar(i1,i0)/TMath::Sqrt(covar(i0,i0)*covar(i1,i1));
3595       if (TMath::Abs(r0-r1)>kEpsilon){
3596         printf("Error 1 - non symetric matrix %d\t%d\t%f",i0,i1,r1-r0);
3597         nerrors++;
3598       }
3599       if (TMath::Abs(r0)>=1){
3600         printf("Error 2 - Wrong correlation %d\t%d\t%f\n",i0,i1,r0);
3601         nerrors++;
3602       }
3603       if (TMath::Abs(r1)>=1){
3604         printf("Error 3 - Wrong correlation %d\t%d\t%f\n",i0,i1,r1);
3605         nerrors++;
3606       }
3607     }
3608   return nerrors;
3609 }
3610
3611
3612
3613 void AliTPCcalibAlign::MakeReportDy(TFile *output){
3614   //
3615   // Draw histogram of dY
3616   //
3617   Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
3618   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
3619
3620   AliTPCcalibAlign *align = this;
3621   TVectorD vecOOP(36);
3622   TVectorD vecOOM(36);
3623   TVectorD vecOIP(36);
3624   TVectorD vecOIM(36);
3625   TVectorD vecOIS(36);
3626   TVectorD vecSec(36);
3627   TCanvas * cOROCdy = new TCanvas("OROC dy","OROC dy",900,600);
3628   cOROCdy->Divide(6,6);
3629   TCanvas * cIROCdy = new TCanvas("IROC dy","IROC dy",900,600);
3630   cIROCdy->Divide(6,6);
3631   TCanvas * cDy = new TCanvas("Dy","Dy",600,700);
3632   cDy->Divide(1,2);
3633   for (Int_t isec=0;isec<36;isec++){
3634     Bool_t isDraw=kFALSE;
3635     vecSec(isec)=isec;
3636     cOROCdy->cd(isec+1);
3637     Int_t secPlus = (isec%18==17)? isec-17:isec+1;
3638     Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
3639     printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
3640     TH1 * hisOOP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus+36);
3641     TH1 * hisOOM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus+36);
3642     TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec);
3643
3644     if (hisOIS) {
3645       hisOIS = (TH1F*)(hisOIS->Clone());
3646       hisOIS->SetDirectory(0);
3647       hisOIS->Scale(1./(hisOIS->GetMaximum()+1));
3648       hisOIS->SetLineColor(kmicolors[0]);
3649       hisOIS->Draw();
3650       isDraw = kTRUE;
3651       vecOIS(isec)=10*hisOIS->GetMean();
3652     }
3653     if (hisOOP) {
3654       hisOOP = (TH1F*)(hisOOP->Clone());
3655       hisOOP->SetDirectory(0);
3656       hisOOP->Scale(1./(hisOOP->GetMaximum()+1));
3657       hisOOP->SetLineColor(kmicolors[1]);
3658       if (isDraw) hisOOP->Draw("same");
3659       if (!isDraw) {hisOOP->Draw(""); isDraw=kTRUE;}
3660       vecOOP(isec)=10*hisOOP->GetMean();
3661     }
3662     if (hisOOM) {
3663       hisOOM = (TH1F*)(hisOOM->Clone());
3664       hisOOM->SetDirectory(0);
3665       hisOOM->Scale(1/(hisOOM->GetMaximum()+1));
3666       hisOOM->SetLineColor(kmicolors[3]);
3667       if (isDraw) hisOOM->Draw("same");
3668       if (!isDraw) {hisOOM->Draw(""); isDraw=kTRUE;}
3669       vecOOM(isec)=10*hisOOM->GetMean();
3670     }
3671   }
3672   //
3673   //
3674   for (Int_t isec=0;isec<36;isec++){
3675     Bool_t isDraw=kFALSE;
3676     cIROCdy->cd(isec+1);
3677     Int_t secPlus = (isec%18==17)? isec-17:isec+1;
3678     Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
3679     printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
3680     TH1 * hisOIP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus);
3681     TH1 * hisOIM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus);
3682     TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec);
3683     if (hisOIS) {
3684       hisOIS = (TH1F*)(hisOIS->Clone());
3685       hisOIS->SetDirectory(0);
3686       hisOIS->Scale(1./(hisOIS->GetMaximum()+1));
3687       hisOIS->SetLineColor(kmicolors[0]);
3688       hisOIS->Draw();
3689       isDraw = kTRUE;
3690       vecOIS(isec)=10*hisOIS->GetMean();
3691     }
3692     if (hisOIP) {
3693       hisOIP = (TH1F*)(hisOIP->Clone());
3694       hisOIP->SetDirectory(0);
3695       hisOIP->Scale(1./(hisOIP->GetMaximum()+1));
3696       hisOIP->SetLineColor(kmicolors[1]);
3697       if (isDraw) hisOIP->Draw("same");
3698       if (!isDraw) {hisOIP->Draw(""); isDraw=kTRUE;}
3699       hisOIP->Draw("same");
3700       vecOIP(isec)=10*hisOIP->GetMean();
3701     }
3702     if (hisOIM) {
3703       hisOIM = (TH1F*)(hisOIM->Clone());
3704       hisOIM->SetDirectory(0);
3705       hisOIM->Scale(1/(hisOIM->GetMaximum()+1));
3706       hisOIM->SetLineColor(kmicolors[3]);
3707       if (isDraw) hisOIM->Draw("same");
3708       if (!isDraw) {hisOIM->Draw(""); isDraw=kTRUE;}
3709       vecOIM(isec)=10*hisOIM->GetMean();
3710     }
3711   }
3712   TGraph* grOIM = new TGraph(36,vecSec.GetMatrixArray(),vecOIM.GetMatrixArray());
3713   TGraph* grOIP = new TGraph(36,vecSec.GetMatrixArray(),vecOIP.GetMatrixArray());
3714   TGraph* grOIS = new TGraph(36,vecSec.GetMatrixArray(),vecOIS.GetMatrixArray());
3715   TGraph* grOOM = new TGraph(36,vecSec.GetMatrixArray(),vecOOM.GetMatrixArray());
3716   TGraph* grOOP = new TGraph(36,vecSec.GetMatrixArray(),vecOOP.GetMatrixArray());
3717   //
3718   grOIS->SetMarkerStyle(kmimarkers[0]);
3719   grOIP->SetMarkerStyle(kmimarkers[1]);
3720   grOIM->SetMarkerStyle(kmimarkers[3]);
3721   grOOP->SetMarkerStyle(kmimarkers[1]);
3722   grOOM->SetMarkerStyle(kmimarkers[3]);
3723   grOIS->SetMarkerColor(kmicolors[0]);
3724   grOIP->SetMarkerColor(kmicolors[1]);
3725   grOIM->SetMarkerColor(kmicolors[3]);
3726   grOOP->SetMarkerColor(kmicolors[1]);
3727   grOOM->SetMarkerColor(kmicolors[3]);
3728   grOIS->SetLineColor(kmicolors[0]);
3729   grOIP->SetLineColor(kmicolors[1]);
3730   grOIM->SetLineColor(kmicolors[3]);
3731   grOOP->SetLineColor(kmicolors[1]);
3732   grOOM->SetLineColor(kmicolors[3]);
3733   grOIS->SetMaximum(1.5);
3734   grOIS->SetMinimum(-1.5);
3735   grOIS->GetXaxis()->SetTitle("Sector number");
3736   grOIS->GetYaxis()->SetTitle("#Delta_{y} (mm)");
3737
3738   cDy->cd(1);
3739   grOIS->Draw("apl");
3740   grOIM->Draw("pl");
3741   grOIP->Draw("pl");
3742   cDy->cd(2);
3743   grOIS->Draw("apl");
3744   grOOM->Draw("pl");
3745   grOOP->Draw("pl");
3746   cOROCdy->SaveAs("picAlign/OROCOROCdy.eps");
3747   cOROCdy->SaveAs("picAlign/OROCOROCdy.gif");
3748   cOROCdy->SaveAs("picAlign/OROCOROCdy.root");
3749   //
3750   cIROCdy->SaveAs("picAlign/OROCIROCdy.eps");
3751   cIROCdy->SaveAs("picAlign/OROCIROCdy.gif");
3752   cIROCdy->SaveAs("picAlign/OROCIROCdy.root");
3753   //
3754   cDy->SaveAs("picAlign/Sectordy.eps");
3755   cDy->SaveAs("picAlign/Sectordy.gif");
3756   cDy->SaveAs("picAlign/Sectordy.root");
3757   if (output){
3758     output->cd();
3759     cOROCdy->Write("OROCOROCDy");
3760     cIROCdy->Write("OROCIROCDy");
3761     cDy->Write("SectorDy");
3762   }
3763 }
3764
3765 void AliTPCcalibAlign::MakeReportDyPhi(TFile */*output*/){
3766   //
3767   //
3768   //
3769   Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
3770   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
3771
3772   AliTPCcalibAlign *align = this;
3773   TCanvas * cOROCdyPhi = new TCanvas("OROC dyphi","OROC dyphi",900,600);
3774   cOROCdyPhi->Divide(6,6);
3775   for (Int_t isec=0;isec<36;isec++){
3776     cOROCdyPhi->cd(isec+1);
3777     Int_t secPlus = (isec%18==17)? isec-17:isec+1;
3778     Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
3779     //printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
3780     TH2F *htemp=0;
3781     TProfile * profdyphiOOP=0,*profdyphiOOM=0,*profdyphiOOS=0;
3782     htemp = (TH2F*) (align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secPlus+36));
3783     if (htemp) profdyphiOOP= htemp->ProfileX();
3784     htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secMinus+36));
3785     if (htemp) profdyphiOOM= htemp->ProfileX();
3786     htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,isec));
3787     if (htemp) profdyphiOOS= htemp->ProfileX();
3788
3789     if (profdyphiOOS){
3790       profdyphiOOS->SetLineColor(kmicolors[0]);
3791       profdyphiOOS->SetMarkerStyle(kmimarkers[0]);
3792       profdyphiOOS->SetMarkerSize(0.2);
3793       profdyphiOOS->SetMaximum(0.5);
3794       profdyphiOOS->SetMinimum(-0.5);
3795       profdyphiOOS->SetXTitle("tan(#phi)");
3796       profdyphiOOS->SetYTitle("#DeltaY (cm)");
3797     }
3798     if (profdyphiOOP){
3799       profdyphiOOP->SetLineColor(kmicolors[1]);
3800       profdyphiOOP->SetMarkerStyle(kmimarkers[1]);
3801       profdyphiOOP->SetMarkerSize(0.2);
3802     }
3803     if (profdyphiOOM){
3804       profdyphiOOM->SetLineColor(kmicolors[3]);
3805       profdyphiOOM->SetMarkerStyle(kmimarkers[3]);
3806       profdyphiOOM->SetMarkerSize(0.2);
3807     }
3808     if (profdyphiOOS){
3809       profdyphiOOS->Draw();
3810     }else{
3811       if (profdyphiOOM) profdyphiOOM->Draw("");
3812       if (profdyphiOOP) profdyphiOOP->Draw("");
3813     }
3814     if (profdyphiOOM) profdyphiOOM->Draw("same");
3815     if (profdyphiOOP) profdyphiOOP->Draw("same");
3816
3817   }
3818 }
3819
3820
3821 //______________________________________________________________________________
3822 void AliTPCcalibAlign::Streamer(TBuffer &R__b)
3823 {
3824   // Stream an object of class AliTPCcalibAlign.
3825   Bool_t isDebug=AliLog::GetDebugLevel("","AliTPCcalibAlign")>0;
3826   if (isDebug) AliSysInfo::SetVerbose(kTRUE);
3827
3828   if (R__b.IsReading()) {
3829     //    R__b.ReadClassBuffer(AliTPCcalibAlign::Class(),this);
3830     UInt_t R__s, R__c;
3831     Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
3832
3833     AliTPCcalibBase::Streamer(R__b); // Stream the base class
3834
3835     // Read the "big data members" from the same Root directory
3836     // Attention: the gDirectory may NOT correspond to the file of the buffer
3837     if (gDirectory){
3838       //      printf("READING from %s\n",gDirectory->GetPath());
3839       for (Int_t i=0; i<2; ++i){
3840         TString hisName=TString::Format("AliTPCcalibAlign.%s.fClusterDelta_%d",GetName(),i);
3841         if (gDirectory->Get(hisName.Data())){
3842           fClusterDelta[i] = dynamic_cast<THn*>(gDirectory->Get((hisName).Data()));
3843         }
3844         if (!fClusterDelta[i]) AliWarning(Form("fClusterDelta[%d] is not read",i));
3845       }
3846       //
3847       for (Int_t i=0; i<4; ++i){
3848         TString hisName=TString::Format("AliTPCcalibAlign.%s.fTrackletDelta_%d",GetName(),i);
3849         if (gDirectory->Get(hisName.Data())){
3850           fTrackletDelta[i] = dynamic_cast<THnSparse*>(gDirectory->Get((hisName).Data()));
3851         }
3852         if (!fTrackletDelta[i]) AliWarning(Form("fTrackletDelta[%d] is not read",i));
3853       }
3854     }
3855     else {
3856       AliError("gDirectory pointer is null");
3857     }
3858     // If the "big data members"were not in the file, try to read them from the object
3859     // This is suppose to restore the backward compatibility
3860     for (Int_t i=0; i<2; ++i) {
3861       if (!fClusterDelta[i]) R__b >> fClusterDelta[i];
3862     }
3863
3864     for (Int_t i=0; i<4; ++i) {
3865       if (!fTrackletDelta[i]) R__b >> fTrackletDelta[i];
3866     }
3867
3868     // Read all the other data members. This is error prone,
3869     // please make sure this part is updated if you change the data members
3870     fDphiHistArray.Streamer(R__b);
3871     fDthetaHistArray.Streamer(R__b);
3872     fDyHistArray.Streamer(R__b);
3873     fDzHistArray.Streamer(R__b);
3874
3875     fDyPhiHistArray.Streamer(R__b);
3876     fDzThetaHistArray.Streamer(R__b);
3877
3878     fDphiZHistArray.Streamer(R__b);
3879     fDthetaZHistArray.Streamer(R__b);
3880     fDyZHistArray.Streamer(R__b);
3881     fDzZHistArray.Streamer(R__b);
3882
3883     fFitterArray12.Streamer(R__b);
3884     fFitterArray9.Streamer(R__b);
3885     fFitterArray6.Streamer(R__b);
3886
3887     fMatrixArray12.Streamer(R__b);
3888     fMatrixArray9.Streamer(R__b);
3889     fMatrixArray6.Streamer(R__b);
3890
3891     fCombinedMatrixArray6.Streamer(R__b);
3892
3893     R__b.ReadFastArray(fPoints,72*72);
3894
3895     R__b >> fNoField;
3896     R__b >> fXIO;
3897     R__b >> fXmiddle;
3898     R__b >> fXquadrant;
3899
3900     fArraySectorIntParam.Streamer(R__b);
3901     fArraySectorIntCovar.Streamer(R__b);
3902
3903     R__b >> fSectorParamA;
3904     R__b >> fSectorCovarA;
3905     R__b >> fSectorParamC;
3906     R__b >> fSectorCovarC;
3907
3908     R__b >> fUseInnerOuter;
3909     R__b >> fgkMergeEntriesCut; // Should we write the static member?
3910   } else {
3911     // R__b.WriteClassBuffer(AliTPCcalibAlign::Class(),this);
3912     if (isDebug) AliSysInfo::AddStamp("aliengStreamer::Start");
3913
3914     R__b.WriteVersion(AliTPCcalibAlign::IsA());
3915     AliTPCcalibBase::Streamer(R__b); // Stream the base class
3916
3917     // Write the "big data members directly in the file in parallel with the object itself
3918     for (Int_t i=0; i<2; ++i){
3919       if (fClusterDelta[i]) fClusterDelta[i]->Write(TString::Format("AliTPCcalibAlign.%s.fClusterDelta_%d",GetName(),i).Data());
3920     }
3921
3922     if (isDebug) AliSysInfo::AddStamp("alignStreamer::fClusterDelta");
3923
3924     for (Int_t i=0; i<4; ++i){
3925       if (fTrackletDelta[i]) fTrackletDelta[i]->Write(TString::Format("AliTPCcalibAlign.%s.fTrackletDelta_%d",GetName(),i).Data());
3926     }
3927
3928     if (isDebug) AliSysInfo::AddStamp("alignStreamer::fTrackletDelta");
3929
3930     // Write all the other data members. This is error prone,
3931     // please make sure this part is updated if you change the data members
3932     fDphiHistArray.Streamer(R__b);
3933     fDthetaHistArray.Streamer(R__b);
3934     fDyHistArray.Streamer(R__b);
3935     fDzHistArray.Streamer(R__b);
3936
3937     fDyPhiHistArray.Streamer(R__b);
3938     fDzThetaHistArray.Streamer(R__b);
3939
3940     fDphiZHistArray.Streamer(R__b);
3941     fDthetaZHistArray.Streamer(R__b);
3942     fDyZHistArray.Streamer(R__b);
3943     fDzZHistArray.Streamer(R__b);
3944
3945     fFitterArray12.Streamer(R__b);
3946     fFitterArray9.Streamer(R__b);
3947     fFitterArray6.Streamer(R__b);
3948
3949     fMatrixArray12.Streamer(R__b);
3950     fMatrixArray9.Streamer(R__b);
3951     fMatrixArray6.Streamer(R__b);
3952
3953     fCombinedMatrixArray6.Streamer(R__b);
3954
3955     R__b.WriteFastArray(fPoints,72*72);
3956
3957     R__b << fNoField;
3958     R__b << fXIO;
3959     R__b << fXmiddle;
3960     R__b << fXquadrant;
3961
3962     fArraySectorIntParam.Streamer(R__b);
3963     fArraySectorIntCovar.Streamer(R__b);
3964
3965     R__b << fSectorParamA;
3966     R__b << fSectorCovarA;
3967     R__b << fSectorParamC;
3968     R__b << fSectorCovarC;
3969
3970     R__b << fUseInnerOuter;
3971     R__b << fgkMergeEntriesCut; // Should we write the static member?
3972     if (isDebug) AliSysInfo::AddStamp("alignStreamer::DefaultStreamer");
3973   }
3974 }