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