TPC/AliTPCcalibAlign.cxx

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