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