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