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