TPC/AliTPCcalibAlign.cxx

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