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