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