TPC/AliTPCcalibAlign.cxx

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