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 "AliExternalTrackParam.h"
 109 #include "AliTPCTracklet.h"
 110 #include "TH1D.h"
 111 #include "TH2F.h"
 112 #include "TVectorD.h"
 113 #include "TTreeStream.h"
 114 #include "TFile.h"
 115 #include "TTree.h"
 116 #include "TF1.h"
 117 #include "TGraphErrors.h"
 118 #include "AliTPCclusterMI.h"
 119 #include "AliTPCseed.h"
 120 #include "AliTracker.h"
 121 #include "TClonesArray.h"
 122 #include "AliExternalComparison.h"
 123
 124
 125 #include "TTreeStream.h"
 126 #include <iostream>
 127 #include <sstream>
 128 using namespace std;
 129
 130 AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0;
 131 ClassImp(AliTPCcalibAlign)
 132
 133
 134
 135
 136 AliTPCcalibAlign* AliTPCcalibAlign::Instance()
 137 {
 138   //
 139   // Singleton implementation
 140   // Returns an instance of this class, it is created if neccessary
 141   //
 142   if (fgInstance == 0){
 143     fgInstance = new AliTPCcalibAlign();
 144   }
 145   return fgInstance;
 146 }
 147
 148
 149
 150
 151 AliTPCcalibAlign::AliTPCcalibAlign()
 152   :  AliTPCcalibBase(),
 153      fDphiHistArray(72*72),
 154      fDthetaHistArray(72*72),
 155      fDyHistArray(72*72),
 156      fDzHistArray(72*72),
 157      //
 158      fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
 159      fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
 160      fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
 161      fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
 162      fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
 163      fDzZHistArray(72*72),        // array of residual histograms  z     -kZz
 164      fFitterArray12(72*72),
 165      fFitterArray9(72*72),
 166      fFitterArray6(72*72),
 167      fMatrixArray12(72*72),
 168      fMatrixArray9(72*72),
 169      fMatrixArray6(72*72),
 170      fCombinedMatrixArray6(72),
 171      fCompTracklet(0),             // tracklet comparison
 172      fNoField(kFALSE)
 173 {
 174   //
 175   // Constructor
 176   //
 177   for (Int_t i=0;i<72*72;++i) {
 178     fPoints[i]=0;
 179   }
 180 }
 181
 182 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
 183   :AliTPCcalibBase(),
 184    fDphiHistArray(72*72),
 185    fDthetaHistArray(72*72),
 186    fDyHistArray(72*72),
 187    fDzHistArray(72*72),
 188    fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
 189    fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
 190    fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
 191    fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
 192    fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
 193    fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     //
 194    fFitterArray12(72*72),
 195    fFitterArray9(72*72),
 196    fFitterArray6(72*72),
 197    fMatrixArray12(72*72),
 198    fMatrixArray9(72*72),
 199    fMatrixArray6(72*72),
 200    fCombinedMatrixArray6(72),
 201    fCompTracklet(0),             // tracklet comparison
 202    fNoField(kFALSE)
 203 {
 204   //
 205   // Constructor
 206   //
 207   SetName(name);
 208   SetTitle(title);
 209   for (Int_t i=0;i<72*72;++i) {
 210     fPoints[i]=0;
 211   }
 212 }
 213
 214
 215 AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
 216   :AliTPCcalibBase(align),
 217    fDphiHistArray(align.fDphiHistArray),
 218    fDthetaHistArray(align.fDthetaHistArray),
 219    fDyHistArray(align.fDyHistArray),
 220    fDzHistArray(align.fDzHistArray),
 221    fDyPhiHistArray(align.fDyPhiHistArray),      // array of residual histograms  y     -kYPhi
 222    fDzThetaHistArray(align.fDzThetaHistArray),    // array of residual histograms  z-z   -kZTheta
 223    fDphiZHistArray(align.fDphiZHistArray),      // array of residual histograms  phi   -kPhiz
 224    fDthetaZHistArray(align.fDthetaZHistArray),    // array of residual histograms  theta -kThetaz
 225    fDyZHistArray(align.fDyZHistArray),        // array of residual histograms  y     -kYz
 226    fDzZHistArray(align.fDzZHistArray),        // array of residual histograms  z     -kZz
 227    //
 228    fFitterArray12(align.fFitterArray12),
 229    fFitterArray9(align.fFitterArray9),
 230    fFitterArray6(align.fFitterArray6),
 231
 232    fMatrixArray12(align.fMatrixArray12),
 233    fMatrixArray9(align.fMatrixArray9),
 234    fMatrixArray6(align.fMatrixArray6),
 235    fCombinedMatrixArray6(align.fCombinedMatrixArray6),
 236    fCompTracklet(align.fCompTracklet),             // tracklet comparison
 237    fNoField(align.fNoField)
 238 {
 239   //
 240   // copy constructor - copy also the content
 241   //
 242   TH1 * his = 0;
 243   TObjArray * arr0=0;
 244   const TObjArray *arr1=0;
 245   for (Int_t index =0; index<72*72; index++){
 246     for (Int_t iarray=0;iarray<10; iarray++){
 247       if (iarray==kY){
 248         arr0 = &fDyHistArray;
 249         arr1 = &align.fDyHistArray;
 250       }
 251       if (iarray==kZ){
 252         arr0 = &fDzHistArray;
 253         arr1 = &align.fDzHistArray;
 254       }
 255       if (iarray==kPhi){
 256         arr0 = &fDphiHistArray;
 257         arr1 = &align.fDphiHistArray;
 258       }
 259       if (iarray==kTheta){
 260         arr0 = &fDthetaHistArray;
 261         arr1 = &align.fDthetaHistArray;
 262       }
 263       if (iarray==kYz){
 264         arr0 = &fDyZHistArray;
 265         arr1 = &align.fDyZHistArray;
 266       }
 267       if (iarray==kZz){
 268         arr0 = &fDzZHistArray;
 269         arr1 = &align.fDzZHistArray;
 270       }
 271       if (iarray==kPhiZ){
 272         arr0 = &fDphiZHistArray;
 273         arr1 = &align.fDphiZHistArray;
 274       }
 275       if (iarray==kThetaZ){
 276         arr0 = &fDthetaZHistArray;
 277         arr1 = &align.fDthetaZHistArray;
 278       }
 279
 280       if (iarray==kYPhi){
 281         arr0 = &fDyPhiHistArray;
 282         arr1 = &align.fDyPhiHistArray;
 283       }
 284       if (iarray==kZTheta){
 285         arr0 = &fDzThetaHistArray;
 286         arr1 = &align.fDzThetaHistArray;
 287       }
 288
 289       if (arr1->At(index)) {
 290         his = (TH1*)arr1->At(index)->Clone();
 291         his->SetDirectory(0);
 292         arr0->AddAt(his,index);
 293       }
 294     }
 295   }
 296 }
 297
 298
 299 AliTPCcalibAlign::~AliTPCcalibAlign() {
 300   //
 301   // destructor
 302   //
 303 }
 304
 305 void AliTPCcalibAlign::Process(AliTPCseed *seed) {
 306   //
 307   //
 308   //
 309
 310   TObjArray tracklets=
 311     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
 312                                     kFALSE,20,2);
 313   //  TObjArray trackletsL=
 314 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
 315 //                                  kFALSE,20,2);
 316 //   TObjArray trackletsQ=
 317 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
 318 //                                  kFALSE,20,2);
 319 //   TObjArray trackletsR=
 320 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
 321 //                                  kFALSE,20,2);
 322   tracklets.SetOwner();
 323  //  trackletsL.SetOwner();
 324 //   trackletsQ.SetOwner();
 325 //   trackletsR.SetOwner();
 326   if (tracklets.GetEntries()==2) {
 327     AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
 328     AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
 329     if (t1->GetSector()>t2->GetSector()) {
 330       AliTPCTracklet* tmp=t1;
 331       t1=t2;
 332       t2=tmp;
 333     }
 334     AliExternalTrackParam *common1=0,*common2=0;
 335     if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
 336       ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
 337     delete common1;
 338     delete common2;
 339   }
 340
 341 }
 342
 343 void AliTPCcalibAlign::Analyze(){
 344   //
 345   // Analyze function
 346   //
 347   EvalFitters();
 348 }
 349
 350
 351 void AliTPCcalibAlign::Terminate(){
 352   //
 353   // Terminate function
 354   // call base terminate + Eval of fitters
 355   //
 356   Info("AliTPCcalibAlign","Terminate");
 357   EvalFitters();
 358   AliTPCcalibBase::Terminate();
 359 }
 360
 361
 362
 363
 364 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
 365                                         const AliExternalTrackParam &tp2,
 366                                         const AliTPCseed * seed,
 367                                         Int_t s1,Int_t s2) {
 368
 369   //
 370   //
 371   //
 372   //
 373   //
 374   // Process function to fill fitters
 375   //
 376   Double_t t1[10],t2[10];
 377   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 378   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 379   x1   =tp1.GetX();
 380   y1   =tp1.GetY();
 381   z1   =tp1.GetZ();
 382   Double_t snp1=tp1.GetSnp();
 383   dydx1=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
 384   Double_t tgl1=tp1.GetTgl();
 385   // dz/dx = 1/(cos(theta)*cos(phi))
 386   dzdx1=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
 387   x2   =tp2.GetX();
 388   y2   =tp2.GetY();
 389   z2   =tp2.GetZ();
 390   Double_t snp2=tp2.GetSnp();
 391   dydx2=snp2/TMath::Sqrt((1.-snp2)*(1.+snp2));
 392   Double_t tgl2=tp2.GetTgl();
 393   dzdx2=tgl2/TMath::Sqrt((1.-snp2)*(1.+snp2));
 394   Int_t accept =   AcceptTracklet(tp1,tp2);
 395   //
 396   //
 397   //
 398   if (fStreamLevel>1 && seed){
 399     TTreeSRedirector *cstream = GetDebugStreamer();
 400     if (cstream){
 401       static TVectorD vec1(5);
 402       static TVectorD vec2(5);
 403       vec1.SetElements(t1);
 404       vec2.SetElements(t2);
 405       AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
 406       AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
 407       (*cstream)<<"Tracklet"<<
 408         "accept="<<accept<<
 409         "run="<<fRun<<              //  run number
 410         "event="<<fEvent<<          //  event number
 411         "time="<<fTime<<            //  time stamp of event
 412         "trigger="<<fTrigger<<      //  trigger
 413         "mag="<<fMagF<<             //  magnetic field
 414         "isOK="<<accept<<           //  flag - used for alignment
 415         "tp1.="<<p1<<
 416         "tp2.="<<p2<<
 417         "v1.="<<&vec1<<
 418         "v2.="<<&vec2<<
 419         "s1="<<s1<<
 420         "s2="<<s2<<
 421         "\n";
 422     }
 423   }
 424   if (accept>0) return;
 425   t1[0]-=134.;
 426   t2[0]-=134.;
 427   t1[5]=0; t2[5]=0;
 428   t1[6]=TMath::Sqrt(tp1.GetSigmaY2()+tp2.GetSigmaY2()); t2[6]=t1[6];
 429   t1[7]=TMath::Sqrt(tp1.GetSigmaZ2()+tp2.GetSigmaZ2()); t2[7]=t1[7];
 430   t1[8]=TMath::Sqrt(tp1.GetSigmaSnp2()+tp2.GetSigmaSnp2()); t2[8]=t1[8];
 431   t1[9]=TMath::Sqrt(tp1.GetSigmaTgl2()+tp2.GetSigmaTgl2()); t2[9]=t1[9];
 432
 433   if (GetDebugLevel()>50) printf("Process track\n");
 434   if (GetDebugLevel()>50) printf("Filling track\n");
 435   //
 436   // fill resolution histograms - previous cut included
 437   if (seed) ProcessDiff(tp1,tp2, seed,s1,s2);
 438   FillHisto(tp1,tp2,s1,s2);
 439   ProcessAlign(t1,t2,s1,s2);
 440 }
 441
 442 void AliTPCcalibAlign::ProcessAlign(Double_t * t1,
 443                                     Double_t * t2,
 444                                     Int_t s1,Int_t s2){
 445   //
 446   // Do intersector alignment
 447   //
 448   Process12(t1,t2,GetOrMakeFitter12(s1,s2));
 449   Process9(t1,t2,GetOrMakeFitter9(s1,s2));
 450   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
 451   ++fPoints[GetIndex(s1,s2)];
 452 }
 453
 454 void AliTPCcalibAlign::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){
 455   //
 456   // Process the debug streamer tree
 457   // Possible to modify selection criteria
 458   // Used with entry list
 459   //
 460   TTreeSRedirector * cstream = new TTreeSRedirector("aligndump.root");
 461
 462   AliTPCcalibAlign *align = this;
 463   //
 464   TVectorD * vec1 = 0;
 465   TVectorD * vec2 = 0;
 466   AliExternalTrackParam * tp1 = 0;
 467   AliExternalTrackParam * tp2 = 0;
 468   Int_t      s1 = 0;
 469   Int_t      s2 = 0;
 470   Int_t npoints =0;
 471   {
 472     Int_t entries=chainTracklet->GetEntries();
 473     for (Int_t i=0; i< entries; i++){
 474       chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
 475       chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
 476       chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
 477       chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
 478       chainTracklet->GetBranch("s1")->SetAddress(&s1);
 479       chainTracklet->GetBranch("s2")->SetAddress(&s2);
 480       chainTracklet->GetEntry(i);
 481       if (!vec1) continue;
 482       if (!vec2) continue;
 483       if (!tp1) continue;
 484       if (!tp2) continue;
 485       if (!vec1->GetMatrixArray()) continue;
 486       if (!vec2->GetMatrixArray()) continue;
 487       // make a local copy
 488       AliExternalTrackParam par1(*tp1);
 489       AliExternalTrackParam par2(*tp2);
 490       TVectorD svec1(*vec1);
 491       TVectorD svec2(*vec2);
 492       //
 493       if (s1==s2) continue;
 494       if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i, npoints,s1,s2);
 495       AliExternalTrackParam  cpar1(par1);
 496       AliExternalTrackParam  cpar2(par2);
 497       Constrain1Pt(cpar1,par2,fNoField);
 498       Constrain1Pt(cpar2,par1,fNoField);
 499       Bool_t acceptComp = kFALSE;
 500       if (comp) acceptComp=comp->AcceptPair(&par1,&par2);
 501       if (comp) acceptComp&=comp->AcceptPair(&cpar1,&cpar2);
 502       //
 503       Int_t reject =   align->AcceptTracklet(par1,par2);
 504       Int_t rejectC =align->AcceptTracklet(cpar1,cpar2);
 505
 506       if (1||fStreamLevel>0){
 507         (*cstream)<<"Tracklet"<<
 508           "s1="<<s1<<
 509           "s2="<<s2<<
 510           "reject="<<reject<<
 511           "rejectC="<<rejectC<<
 512           "acceptComp="<<acceptComp<<
 513           "tp1.="<<&par1<<
 514           "tp2.="<<&par2<<
 515           "ctp1.="<<&cpar1<<
 516           "ctp2.="<<&cpar2<<
 517           "v1.="<<&svec1<<
 518           "v2.="<<&svec2<<
 519           "\n";
 520       }
 521       //
 522       if (fNoField){
 523         //
 524         //
 525       }
 526       if (acceptComp) comp->Process(&cpar1,&cpar2);
 527       //
 528       if (reject>0 || rejectC>0) continue;
 529       npoints++;
 530       align->ProcessTracklets(cpar1,cpar2,0,s1,s2);
 531       align->ProcessTracklets(cpar2,cpar1,0,s2,s1);
 532     }
 533   }
 534   delete cstream;
 535 }
 536
 537
 538 Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
 539                                        const AliExternalTrackParam &p2){
 540
 541   //
 542   // Accept pair of tracklets?
 543   //
 544   /*
 545   // resolution cuts
 546   TCut cutS0("sqrt(tp2.fC[0]+tp1.fC[0])<0.2");
 547   TCut cutS1("sqrt(tp2.fC[2]+tp1.fC[2])<0.2");
 548   TCut cutS2("sqrt(tp2.fC[5]+tp1.fC[5])<0.01");
 549   TCut cutS3("sqrt(tp2.fC[9]+tp1.fC[9])<0.01");
 550   TCut cutS4("sqrt(tp2.fC[14]+tp1.fC[14])<0.25");
 551   TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
 552   //
 553   // parameters matching cuts
 554   TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
 555   TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
 556   TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03");
 557   TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03");
 558   TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
 559   TCut cutPP4("abs(tp1.fP[4]-tp2.fP[4])/sqrt(tp2.fC[14]+tp1.fC[14])<3");
 560   TCut cutP=cutP0+cutP1+cutP2+cutP3+cutP4+cutPP4;
 561   */
 562   //
 563   // resolution cuts
 564   Int_t reject=0;
 565   const Double_t *cp1 = p1.GetCovariance();
 566   const Double_t *cp2 = p2.GetCovariance();
 567   if (TMath::Sqrt(cp1[0]+cp2[0])>0.2)  reject|=1;;
 568   if (TMath::Sqrt(cp1[2]+cp2[2])>0.2)  reject|=2;
 569   if (TMath::Sqrt(cp1[5]+cp2[5])>0.01) reject|=4;
 570   if (TMath::Sqrt(cp1[9]+cp2[9])>0.01) reject|=8;
 571   if (TMath::Sqrt(cp1[14]+cp2[14])>0.2) reject|=16;
 572
 573   //parameters difference
 574   const Double_t *tp1 = p1.GetParameter();
 575   const Double_t *tp2 = p2.GetParameter();
 576   if (TMath::Abs(tp1[0]-tp2[0])>0.6) reject|=32;
 577   if (TMath::Abs(tp1[1]-tp2[1])>0.6) reject|=64;
 578   if (TMath::Abs(tp1[2]-tp2[2])>0.03) reject|=128;
 579   if (TMath::Abs(tp1[3]-tp2[3])>0.03) reject|=526;
 580   if (TMath::Abs(tp1[4]-tp2[4])>0.4) reject|=1024;
 581   if (TMath::Abs(tp1[4]-tp2[4])/TMath::Sqrt(cp1[14]+cp2[14])>4) reject|=2048;
 582
 583   //
 584   if (TMath::Abs(tp2[1])>235) reject|=2*4096;
 585
 586   if (fNoField){
 587
 588   }
 589
 590   return reject;
 591 }
 592
 593
 594
 595 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
 596                                     const AliExternalTrackParam &t2,
 597                                     const AliTPCseed *seed,
 598                                     Int_t s1,Int_t s2)
 599 {
 600   //
 601   // Process local residuals function
 602   //
 603   TVectorD vecX(160);
 604   TVectorD vecY(160);
 605   TVectorD vecZ(160);
 606   TVectorD vecClY(160);
 607   TVectorD vecClZ(160);
 608   TClonesArray arrCl("AliTPCclusterMI",160);
 609   arrCl.ExpandCreateFast(160);
 610   Int_t count1=0, count2=0;
 611
 612   for (Int_t i=0;i<160;++i) {
 613     AliTPCclusterMI *c=seed->GetClusterPointer(i);
 614     vecX[i]=0;
 615     vecY[i]=0;
 616     vecZ[i]=0;
 617     if (!c) continue;
 618     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
 619     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
 620     vecClY[i] = c->GetY();
 621     vecClZ[i] = c->GetZ();
 622     cl=*c;
 623     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
 624     if (c->GetDetector()==s1) ++count1;
 625     if (c->GetDetector()==s2) ++count2;
 626     Double_t gxyz[3],xyz[3];
 627     t1.GetXYZ(gxyz);
 628     Float_t bz = AliTracker::GetBz(gxyz);
 629     par->GetYAt(c->GetX(), bz, xyz[1]);
 630     par->GetZAt(c->GetX(), bz, xyz[2]);
 631     vecX[i] = c->GetX();
 632     vecY[i]= xyz[1];
 633     vecZ[i]= xyz[2];
 634   }
 635   //
 636   //
 637   if (fStreamLevel>5){
 638     //
 639     // huge output - cluster residuals to be investigated
 640     //
 641     TTreeSRedirector *cstream = GetDebugStreamer();
 642     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
 643     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
 644     /*
 645
 646       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");
 647
 648     */
 649
 650     if (cstream){
 651       (*cstream)<<"Track"<<
 652         "run="<<fRun<<              //  run number
 653         "event="<<fEvent<<          //  event number
 654         "time="<<fTime<<            //  time stamp of event
 655         "trigger="<<fTrigger<<      //  trigger
 656         "mag="<<fMagF<<             //  magnetic field
 657         "Cl.="<<&arrCl<<
 658         //"tp0.="<<p0<<
 659         "tp1.="<<p1<<
 660         "tp2.="<<p2<<
 661         "vtX.="<<&vecX<<
 662         "vtY.="<<&vecY<<
 663         "vtZ.="<<&vecZ<<
 664         "vcY.="<<&vecClY<<
 665         "vcZ.="<<&vecClZ<<
 666         "s1="<<s1<<
 667         "s2="<<s2<<
 668         "c1="<<count1<<
 669         "c2="<<count2<<
 670         "\n";
 671     }
 672   }
 673 }
 674
 675
 676
 677
 678 void AliTPCcalibAlign::Process12(const Double_t *t1,
 679                                  const Double_t *t2,
 680                                  TLinearFitter *fitter) {
 681   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 682   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 683   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 684   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 685   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 686   //
 687   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
 688   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
 689   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
 690
 691
 692
 693   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 694   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 695
 696   //
 697   Double_t sy    = t1[6];
 698   Double_t sz    = t1[7];
 699   Double_t sdydx = t1[8];
 700   Double_t sdzdx = t1[9];
 701
 702   Double_t p[12];
 703   Double_t value;
 704
 705   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 706   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 707   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 708   for (Int_t i=0; i<12;i++) p[i]=0.;
 709   p[3+0] = x1;          // a10
 710   p[3+1] = y1;          // a11
 711   p[3+2] = z1;          // a12
 712   p[9+1] = 1.;          // a13
 713   p[0+1] = y1*dydx2;    // a01
 714   p[0+2] = z1*dydx2;    // a02
 715   p[9+0] = dydx2;       // a03
 716   value  = y2;
 717   fitter->AddPoint(p,value,sy);
 718
 719   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 720   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 721   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 722   for (Int_t i=0; i<12;i++) p[i]=0.;
 723   p[6+0] = x1;           // a20
 724   p[6+1] = y1;           // a21
 725   p[6+2] = z1;           // a22
 726   p[9+2] = 1.;           // a23
 727   p[0+1] = y1*dzdx2;     // a01
 728   p[0+2] = z1*dzdx2;     // a02
 729   p[9+0] = dzdx2;        // a03
 730   value  = z2;
 731   fitter->AddPoint(p,value,sz);
 732
 733   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 734   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 735   for (Int_t i=0; i<12;i++) p[i]=0.;
 736   p[3+0] = 1.;           // a10
 737   p[3+1] = dydx1;        // a11
 738   p[3+2] = dzdx1;        // a12
 739   p[0+0] = -dydx2;       // a00
 740   p[0+1] = -dydx1*dydx2; // a01
 741   p[0+2] = -dzdx1*dydx2; // a02
 742   value  = 0.;
 743   fitter->AddPoint(p,value,sdydx);
 744
 745   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 746   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 747   for (Int_t i=0; i<12;i++) p[i]=0.;
 748   p[6+0] = 1;            // a20
 749   p[6+1] = dydx1;        // a21
 750   p[6+2] = dzdx1;        // a22
 751   p[0+0] = -dzdx2;       // a00
 752   p[0+1] = -dydx1*dzdx2; // a01
 753   p[0+2] = -dzdx1*dzdx2; // a02
 754   value  = 0.;
 755   fitter->AddPoint(p,value,sdzdx);
 756 }
 757
 758 void AliTPCcalibAlign::Process9(Double_t *t1,
 759                                 Double_t *t2,
 760                                 TLinearFitter *fitter) {
 761   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 762   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 763   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 764   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 765   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 766   //
 767   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
 768   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
 769   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8]
 770
 771
 772   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 773   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 774   //
 775   Double_t sy    = t1[6];
 776   Double_t sz    = t1[7];
 777   Double_t sdydx = t1[8];
 778   Double_t sdzdx = t1[9];
 779   //
 780   Double_t p[12];
 781   Double_t value;
 782
 783   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 784   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 785   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 786   for (Int_t i=0; i<12;i++) p[i]=0.;
 787   p[2]   += x1;           // a10
 788   //p[]  +=1;             // a11
 789   p[3]   += z1;           // a12
 790   p[7]   += 1;            // a13
 791   p[0]   += y1*dydx2;     // a01
 792   p[1]   += z1*dydx2;     // a02
 793   p[6]   += dydx2;        // a03
 794   value   = y2-y1;        //-a11
 795   fitter->AddPoint(p,value,sy);
 796   //
 797   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 798   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 799   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 800   for (Int_t i=0; i<12;i++) p[i]=0.;
 801   p[4]   += x1;           // a20
 802   p[5]   += y1;           // a21
 803   //p[]  += z1;           // a22
 804   p[8]   += 1.;           // a23
 805   p[0]   += y1*dzdx2;     // a01
 806   p[1]   += z1*dzdx2;     // a02
 807   p[6]   += dzdx2;        // a03
 808   value  = z2-z1;         //-a22
 809   fitter->AddPoint(p,value,sz);
 810
 811   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 812   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 813   for (Int_t i=0; i<12;i++) p[i]=0.;
 814   p[2]   += 1.;           // a10
 815   //p[]  += dydx1;      // a11
 816   p[3]   += dzdx1;        // a12
 817   //p[]  += -dydx2;       // a00
 818   p[0]   += -dydx1*dydx2; // a01
 819   p[1]   += -dzdx1*dydx2; // a02
 820   value  = -dydx1+dydx2;  // -a11 + a00
 821   fitter->AddPoint(p,value,sdydx);
 822
 823   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 824   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 825   for (Int_t i=0; i<12;i++) p[i]=0.;
 826   p[4]   += 1;            // a20
 827   p[5]   += dydx1;        // a21
 828   //p[]  += dzdx1;        // a22
 829   //p[]  += -dzdx2;       // a00
 830   p[0]   += -dydx1*dzdx2; // a01
 831   p[1]   += -dzdx1*dzdx2; // a02
 832   value  = -dzdx1+dzdx2;  // -a22 + a00
 833   fitter->AddPoint(p,value,sdzdx);
 834 }
 835
 836 void AliTPCcalibAlign::Process6(Double_t *t1,
 837                                 Double_t *t2,
 838                                 TLinearFitter *fitter) {
 839   // x2    =  1  *x1 +-a01*y1 + 0      +a03
 840   // y2    =  a01*x1 + 1  *y1 + 0      +a13
 841   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
 842   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 843   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 844   //
 845   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
 846   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
 847   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5]
 848
 849   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 850   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 851
 852   //
 853   Double_t sy    = t1[6];
 854   Double_t sz    = t1[7];
 855   Double_t sdydx = t1[8];
 856   Double_t sdzdx = t1[9];
 857
 858   Double_t p[12];
 859   Double_t value;
 860   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 861   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 862   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 863   for (Int_t i=0; i<12;i++) p[i]=0.;
 864   p[0]   += x1;           // a10
 865   //p[]  +=1;             // a11
 866   //p[]  += z1;           // a12
 867   p[4]   += 1;            // a13
 868   p[0]   += -y1*dydx2;    // a01
 869   //p[]  += z1*dydx2;     // a02
 870   p[3]   += dydx2;        // a03
 871   value   = y2-y1;        //-a11
 872   fitter->AddPoint(p,value,sy);
 873   //
 874   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 875   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 876   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 877   for (Int_t i=0; i<12;i++) p[i]=0.;
 878   p[1]   += x1;           // a20
 879   p[2]   += y1;           // a21
 880   //p[]  += z1;           // a22
 881   p[5]   += 1.;           // a23
 882   p[0]   += -y1*dzdx2;    // a01
 883   //p[]   += z1*dzdx2;     // a02
 884   p[3]   += dzdx2;        // a03
 885   value  = z2-z1;         //-a22
 886   fitter->AddPoint(p,value,sz);
 887
 888   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 889   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 890   for (Int_t i=0; i<12;i++) p[i]=0.;
 891   p[0]   += 1.;           // a10
 892   //p[]  += dydx1;      // a11
 893   //p[]   += dzdx1;        // a12
 894   //p[]  += -dydx2;       // a00
 895   //p[0]   +=  dydx1*dydx2; // a01         FIXME- 0912 MI
 896   //p[]   += -dzdx1*dydx2; // a02
 897   value  = -dydx1+dydx2;  // -a11 + a00
 898   fitter->AddPoint(p,value,sdydx);
 899
 900   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 901   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 902   for (Int_t i=0; i<12;i++) p[i]=0.;
 903   p[1]   += 1;            // a20
 904   //  p[2]   += dydx1;        // a21   FIXME- 0912 MI
 905   //p[]  += dzdx1;        // a22
 906   //p[]  += -dzdx2;       // a00
 907   //p[0]   +=  dydx1*dzdx2; // a01     FIXME- 0912 MI
 908   //p[]  += -dzdx1*dzdx2; // a02
 909   value  = -dzdx1+dzdx2;  // -a22 + a00
 910   fitter->AddPoint(p,value,sdzdx);
 911 }
 912
 913
 914
 915
 916 void AliTPCcalibAlign::EvalFitters() {
 917   //
 918   // Analyze function
 919   //
 920   // Perform the fitting using linear fitters
 921   //
 922   Int_t kMinPoints =50;
 923   TLinearFitter *f;
 924   TFile fff("alignDebug.root","recreate");
 925   for (Int_t s1=0;s1<72;++s1)
 926     for (Int_t s2=0;s2<72;++s2){
 927       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 928         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 929         if (f->Eval()!=0) {
 930           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 931           f->Write(Form("f12_%d_%d",s1,s2));
 932         }else{
 933           f->Write(Form("f12_%d_%d",s1,s2));
 934         }
 935       }
 936       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 937         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 938         if (f->Eval()!=0) {
 939           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 940         }else{
 941           f->Write(Form("f9_%d_%d",s1,s2));
 942         }
 943       }
 944       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 945         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 946         if (f->Eval()!=0) {
 947           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 948         }else{
 949           f->Write(Form("f6_%d_%d",s1,s2));
 950         }
 951       }
 952     }
 953   TMatrixD mat(4,4);
 954   for (Int_t s1=0;s1<72;++s1)
 955     for (Int_t s2=0;s2<72;++s2){
 956       if (GetTransformation12(s1,s2,mat)){
 957         fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2));
 958       }
 959       if (GetTransformation9(s1,s2,mat)){
 960         fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2));
 961       }
 962       if (GetTransformation6(s1,s2,mat)){
 963         fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2));
 964       }
 965     }
 966   //this->Write("align");
 967
 968 }
 969
 970 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
 971   //
 972   // get or make fitter - general linear transformation
 973   //
 974   static Int_t counter12=0;
 975   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]");
 976   TLinearFitter * fitter = GetFitter12(s1,s2);
 977   if (fitter) return fitter;
 978   //  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]");
 979   fitter =new TLinearFitter(&f12,"");
 980   fitter->StoreData(kFALSE);
 981   fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
 982   counter12++;
 983   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
 984   return fitter;
 985 }
 986
 987 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
 988   //
 989   //get or make fitter - general linear transformation - no scaling
 990   //
 991   static Int_t counter9=0;
 992   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 993   TLinearFitter * fitter = GetFitter9(s1,s2);
 994   if (fitter) return fitter;
 995   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 996   fitter =new TLinearFitter(&f9,"");
 997   fitter->StoreData(kFALSE);
 998   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
 999   counter9++;
1000   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
1001   return fitter;
1002 }
1003
1004 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
1005   //
1006   // get or make fitter  - 6 paramater linear tranformation
1007   //                     - no scaling
1008   //                     - rotation x-y
1009   //                     - tilting x-z, y-z
1010   static Int_t counter6=0;
1011   static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1012   TLinearFitter * fitter = GetFitter6(s1,s2);
1013   if (fitter) return fitter;
1014   //  fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1015   fitter=new TLinearFitter(&f6,"");
1016   fitter->StoreData(kTRUE);
1017   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
1018   counter6++;
1019   if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
1020   return fitter;
1021 }
1022
1023
1024
1025
1026
1027 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
1028   //
1029   // GetTransformation matrix - 12 paramaters - generael linear transformation
1030   //
1031   if (!GetFitter12(s1,s2))
1032     return false;
1033   else {
1034     TVectorD p(12);
1035     GetFitter12(s1,s2)->GetParameters(p);
1036     a.ResizeTo(4,4);
1037     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
1038     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
1039     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
1040     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1041     return true;
1042   }
1043 }
1044
1045 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
1046   //
1047   // GetTransformation matrix - 9 paramaters - general linear transformation
1048   //                            No scaling
1049   //
1050   if (!GetFitter9(s1,s2))
1051     return false;
1052   else {
1053     TVectorD p(9);
1054     GetFitter9(s1,s2)->GetParameters(p);
1055     a.ResizeTo(4,4);
1056     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
1057     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
1058     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
1059     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1060     return true;
1061   }
1062 }
1063
1064 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
1065   //
1066   // GetTransformation matrix - 6  paramaters
1067   //                            3  translation
1068   //                            1  rotation -x-y
1069   //                            2  tilting x-z y-z
1070   if (!GetFitter6(s1,s2))
1071     return false;
1072   else {
1073     TVectorD p(6);
1074     GetFitter6(s1,s2)->GetParameters(p);
1075     a.ResizeTo(4,4);
1076     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
1077     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
1078     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
1079     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1080     return true;
1081   }
1082 }
1083
1084 void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
1085                                         const AliExternalTrackParam &tp2,
1086                                         Int_t s1,Int_t s2) {
1087   //
1088   // Fill residual histograms
1089   // Innner-Outer
1090   // Left right - x-y
1091   // A-C side
1092   if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0)  {
1093     GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
1094     GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1095     GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
1096     GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
1097     //
1098     GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
1099     GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1100     GetHisto(kYz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetY()-tp2.GetY());
1101     GetHisto(kZz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetZ()-tp2.GetZ());
1102     //
1103     GetHisto(kYPhi,s1,s2,kTRUE)->Fill(tp1.GetSnp(),tp1.GetY()-tp2.GetY());
1104     GetHisto(kZTheta,s1,s2,kTRUE)->Fill(tp1.GetTgl(),tp1.GetZ()-tp2.GetZ());
1105
1106
1107   }
1108 }
1109
1110
1111
1112 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
1113 {
1114   //
1115   // return specified residual histogram - it is only QA
1116   // if force specified the histogram and given histogram is not existing
1117   //  new histogram is created
1118   //
1119   if (GetIndex(s1,s2)>=72*72) return 0;
1120   TObjArray *histoArray=0;
1121   switch (type) {
1122   case kY:
1123     histoArray = &fDyHistArray; break;
1124   case kZ:
1125     histoArray = &fDzHistArray; break;
1126   case kPhi:
1127     histoArray = &fDphiHistArray; break;
1128   case kTheta:
1129     histoArray = &fDthetaHistArray; break;
1130   case kYPhi:
1131     histoArray = &fDyPhiHistArray; break;
1132   case kZTheta:
1133     histoArray = &fDzThetaHistArray; break;
1134   case kYz:
1135     histoArray = &fDyZHistArray; break;
1136   case kZz:
1137     histoArray = &fDzZHistArray; break;
1138   case kPhiZ:
1139     histoArray = &fDphiZHistArray; break;
1140   case kThetaZ:
1141     histoArray = &fDthetaZHistArray; break;
1142   }
1143   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
1144   if (histo) return histo;
1145   if (force==kFALSE) return 0;
1146   //
1147   stringstream name;
1148   stringstream title;
1149   switch (type) {
1150   case kY:
1151     name<<"hist_y_"<<s1<<"_"<<s2;
1152     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1153     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1154     break;
1155   case kZ:
1156     name<<"hist_z_"<<s1<<"_"<<s2;
1157     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1158     histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1159     break;
1160   case kPhi:
1161     name<<"hist_phi_"<<s1<<"_"<<s2;
1162     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1163     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1164     break;
1165   case kTheta:
1166     name<<"hist_theta_"<<s1<<"_"<<s2;
1167     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1168     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1169     break;
1170     //
1171     //
1172   case kYPhi:
1173     name<<"hist_yphi_"<<s1<<"_"<<s2;
1174     title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1175     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,128,-0.3,0.3); // +/- 3 mm
1176     break;
1177   case kZTheta:
1178     name<<"hist_ztheta_"<<s1<<"_"<<s2;
1179     title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1180     histo = new TH2F(name.str().c_str(),title.str().c_str(),128,20,-1,1,-0.3,0.3); // +/- 3 mm
1181     break;
1182     //
1183     //
1184     //
1185   case kYz:
1186     name<<"hist_yz_"<<s1<<"_"<<s2;
1187     title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1188     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1189     break;
1190   case kZz:
1191     name<<"hist_zz_"<<s1<<"_"<<s2;
1192     title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1193     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1194     break;
1195   case kPhiZ:
1196     name<<"hist_phiz_"<<s1<<"_"<<s2;
1197     title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1198     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1199     break;
1200   case kThetaZ:
1201     name<<"hist_thetaz_"<<s1<<"_"<<s2;
1202     title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1203     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1204     break;
1205
1206
1207   }
1208   histo->SetDirectory(0);
1209   histoArray->AddAt(histo,GetIndex(s1,s2));
1210   return histo;
1211 }
1212
1213 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
1214                                            Int_t i0, Int_t i1, FitType type)
1215 {
1216   //
1217   //
1218   //
1219   TMatrixD mat;
1220   //TObjArray *fitArray=0;
1221   Double_t xsec[1000];
1222   Double_t ysec[1000];
1223   Int_t npoints=0;
1224   for (Int_t isec = sec0; isec<=sec1; isec++){
1225     Int_t isec2 = (isec+dsec)%72;
1226     switch (type) {
1227     case k6:
1228       GetTransformation6(isec,isec2,mat);break;
1229     case k9:
1230       GetTransformation9(isec,isec2,mat);break;
1231     case k12:
1232       GetTransformation12(isec,isec2,mat);break;
1233     }
1234     xsec[npoints]=isec;
1235     ysec[npoints]=mat(i0,i1);
1236     ++npoints;
1237   }
1238   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1239   Char_t name[1000];
1240   sprintf(name,"Mat[%d,%d]  Type=%d",i0,i1,type);
1241   gr->SetName(name);
1242   return gr;
1243 }
1244
1245 void  AliTPCcalibAlign::MakeTree(const char *fname){
1246   //
1247   // make tree with alignment cosntant  -
1248   // For  QA visualization
1249   //
1250   /*
1251    TFile f("CalibObjects.root");
1252    TObjArray *array  = (TObjArray*)f.Get("TPCCalib");
1253    AliTPCcalibAlign   *alignTPC = (AliTPCcalibAlign   *)array->At(0);
1254    alignTPC->MakeTree("alignTree.root");
1255    TFile falign("alignTree.root");
1256    Align->Draw("dy")
1257    */
1258   const Int_t kMinPoints=50;
1259   TTreeSRedirector cstream(fname);
1260   for (Int_t s1=0;s1<72;++s1)
1261     for (Int_t s2=0;s2<72;++s2){
1262       TMatrixD m6;
1263       TMatrixD m6FX;
1264       TMatrixD m9;
1265       TMatrixD m12;
1266       Double_t dy=0, dz=0, dphi=0,dtheta=0;
1267       Double_t sy=0, sz=0, sphi=0,stheta=0;
1268       Double_t ny=0, nz=0, nphi=0,ntheta=0;
1269       Double_t chi2v12=0, chi2v9=0, chi2v6=0;
1270       Int_t npoints=0;
1271       TLinearFitter * fitter = 0;
1272       if (fPoints[GetIndex(s1,s2)]>kMinPoints){
1273         //
1274         //
1275         //
1276         fitter = GetFitter12(s1,s2);
1277         npoints = fitter->GetNpoints();
1278         chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1279         //
1280         fitter = GetFitter9(s1,s2);
1281         npoints = fitter->GetNpoints();
1282         chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1283         //
1284         fitter = GetFitter6(s1,s2);
1285         npoints = fitter->GetNpoints();
1286         chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1287
1288         //
1289         GetTransformation6(s1,s2,m6);
1290         GetTransformation9(s1,s2,m9);
1291         GetTransformation12(s1,s2,m12);
1292         //
1293         fitter = GetFitter6(s1,s2);
1294         fitter->FixParameter(3,0);
1295         fitter->Eval();
1296         GetTransformation6(s1,s2,m6FX);
1297         //
1298         TH1 * his=0;
1299         his = GetHisto(kY,s1,s2);
1300         if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1301         his = GetHisto(kZ,s1,s2);
1302         if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1303         his = GetHisto(kPhi,s1,s2);
1304         if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
1305         his = GetHisto(kTheta,s1,s2);
1306         if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
1307         //
1308
1309       }
1310
1311       // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1312       // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1313       // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1314       // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1315       // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1316       //
1317       //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
1318       //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
1319       //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
1320
1321       //
1322       //
1323       // dy:-(134*m6.fElements[4]+m6.fElements[7])
1324       //
1325       // dphi:-(m6.fElements[4])
1326       //
1327       // dz:134*m6.fElements[8]+m6.fElements[11]
1328       //
1329       // dtheta:m6.fElements[8]
1330       //
1331       cstream<<"Align"<<
1332         "s1="<<s1<<     // reference sector
1333         "s2="<<s2<<     // sector to align
1334         "m6FX.="<<&m6FX<<   // tranformation matrix
1335         "m6.="<<&m6<<   // tranformation matrix
1336         "m9.="<<&m9<<   //
1337         "m12.="<<&m12<<
1338         "chi2v12="<<chi2v12<<
1339         "chi2v9="<<chi2v9<<
1340         "chi2v6="<<chi2v6<<
1341         //               histograms mean RMS and entries
1342         "dy="<<dy<<
1343         "sy="<<sy<<
1344         "ny="<<ny<<
1345         "dz="<<dz<<
1346         "sz="<<sz<<
1347         "nz="<<nz<<
1348         "dphi="<<dphi<<
1349         "sphi="<<sphi<<
1350         "nphi="<<nphi<<
1351         "dtheta="<<dtheta<<
1352         "stheta="<<stheta<<
1353         "ntheta="<<ntheta<<
1354         "\n";
1355     }
1356
1357 }
1358
1359
1360 //_____________________________________________________________________
1361 Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
1362   //
1363   // merge function
1364   //
1365   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
1366   if (!list)
1367     return 0;
1368   if (list->IsEmpty())
1369     return 1;
1370
1371   TIterator* iter = list->MakeIterator();
1372   TObject* obj = 0;
1373   iter->Reset();
1374   Int_t count=0;
1375   //
1376   TString str1(GetName());
1377   while((obj = iter->Next()) != 0)
1378     {
1379       AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
1380       if (entry == 0) continue;
1381       if (str1.CompareTo(entry->GetName())!=0) continue;
1382       Add(entry);
1383       count++;
1384     }
1385   return count;
1386 }
1387
1388
1389 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
1390   //
1391   // Add entry - used for merging of compoents
1392   //
1393   for (Int_t i=0; i<72;i++){
1394     for (Int_t j=0; j<72;j++){
1395       if (align->fPoints[GetIndex(i,j)]<10) continue;
1396       fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
1397       //
1398       //
1399       //
1400       for (Int_t itype=0; itype<10; itype++){
1401         TH1 * his0=0, *his1=0;
1402         his0 = GetHisto((HistoType)itype,i,j);
1403         his1 = align->GetHisto((HistoType)itype,i,j);
1404         if (his1){
1405           if (his0) his0->Add(his1);
1406           else {
1407             his0 = GetHisto(kY,i,j,kTRUE);
1408             his0->Add(his1);
1409           }
1410         }
1411       }
1412     }
1413   }
1414   TLinearFitter *f0=0;
1415   TLinearFitter *f1=0;
1416   for (Int_t i=0; i<72;i++){
1417     for (Int_t j=0; j<72;j++){
1418       if (align->fPoints[GetIndex(i,j)]<20) continue;
1419       //
1420       //
1421       // fitter12
1422       f0 =  GetFitter12(i,j);
1423       f1 =  align->GetFitter12(i,j);
1424       if (f1 &&f1->Eval()){
1425         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1426         else {
1427           fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
1428         }
1429       }
1430       //
1431       // fitter9
1432       f0 =  GetFitter9(i,j);
1433       f1 =  align->GetFitter9(i,j);
1434       if (f1&&f1->Eval()){
1435         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1436         else {
1437           fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
1438         }
1439       }
1440       f0 =  GetFitter6(i,j);
1441       f1 =  align->GetFitter6(i,j);
1442       if (f1 &&f1->Eval()){
1443         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1444         else {
1445           fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
1446         }
1447       }
1448     }
1449   }
1450 }
1451
1452 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){
1453   //
1454   // GetTransformed value
1455   //
1456   //
1457   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1458   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1459   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1460   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1461   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1462
1463
1464   const TMatrixD * mat = GetTransformation(s1,s2,type);
1465   if (!mat) {
1466     if (value==0) return x1;
1467     if (value==1) return y1;
1468     if (value==2) return z1;
1469     if (value==3) return dydx1;
1470     if (value==4) return dzdx1;
1471     //
1472     if (value==5) return dydx1;
1473     if (value==6) return dzdx1;
1474     return 0;
1475   }
1476   Double_t valT=0;
1477
1478   if (value==0){
1479     valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
1480   }
1481
1482   if (value==1){
1483     valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
1484   }
1485   if (value==2){
1486     valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
1487   }
1488   if (value==3){
1489     //    dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1490     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
1491     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1492   }
1493
1494   if (value==4){
1495     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1496     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
1497     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1498   }
1499   //
1500   if (value==5){
1501     // onlys shift in angle
1502     //    dydx2 =  (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1503     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1;
1504   }
1505
1506   if (value==6){
1507     // only shift in angle
1508     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1509     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1;
1510   }
1511   //
1512   return valT;
1513 }
1514
1515
1516 void  AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){
1517   //
1518   // Update track parameters t1
1519   //
1520   TMatrixD vecXk(5,1);    // X vector
1521   TMatrixD covXk(5,5);    // X covariance
1522   TMatrixD matHk(1,5);    // vector to mesurement
1523   TMatrixD measR(1,1);    // measurement error
1524   //TMatrixD matQk(5,5);    // prediction noise vector
1525   TMatrixD vecZk(1,1);    // measurement
1526   //
1527   TMatrixD vecYk(1,1);    // Innovation or measurement residual
1528   TMatrixD matHkT(5,1);
1529   TMatrixD matSk(1,1);    // Innovation (or residual) covariance
1530   TMatrixD matKk(5,1);    // Optimal Kalman gain
1531   TMatrixD mat1(5,5);     // update covariance matrix
1532   TMatrixD covXk2(5,5);   //
1533   TMatrixD covOut(5,5);
1534   //
1535   Double_t *param1=(Double_t*) track1.GetParameter();
1536   Double_t *covar1=(Double_t*) track1.GetCovariance();
1537
1538   //
1539   // copy data to the matrix
1540   for (Int_t ipar=0; ipar<5; ipar++){
1541     vecXk(ipar,0) = param1[ipar];
1542     for (Int_t jpar=0; jpar<5; jpar++){
1543       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
1544     }
1545   }
1546   //
1547   //
1548   //
1549   vecZk(0,0) = track2.GetParameter()[4];   // 1/pt measurement from track 2
1550   measR(0,0) = track2.GetCovariance()[14];  // 1/pt measurement error
1551   if (noField) {
1552     measR(0,0)*=0.000000001;
1553     vecZk(0,0)=0.;
1554   }
1555   //
1556   matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0;
1557   matHk(0,3)= 0;    matHk(0,4)= 1;           // vector to measurement
1558   //
1559   //
1560   //
1561   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
1562   matHkT=matHk.T(); matHk.T();
1563   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
1564   matSk.Invert();
1565   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
1566   vecXk += matKk*vecYk;                      //  updated vector
1567   mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
1568   covXk2 = (mat1-(matKk*matHk));
1569   covOut =  covXk2*covXk;
1570   //
1571   //
1572   //
1573   // copy from matrix to parameters
1574   if (0) {
1575     covOut.Print();
1576     vecXk.Print();
1577     covXk.Print();
1578     track1.Print();
1579     track2.Print();
1580   }
1581
1582   for (Int_t ipar=0; ipar<5; ipar++){
1583     param1[ipar]= vecXk(ipar,0) ;
1584     for (Int_t jpar=0; jpar<5; jpar++){
1585       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
1586     }
1587   }
1588
1589 }
1590
1591 void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){
1592   //
1593   //  Global Align -combine the partial alignment of pair of sectors
1594   //  minPoints - minimal number of points - don't use sector alignment wit smaller number
1595   //  sysError  - error added to the alignemnt error
1596   //
1597   AliTPCcalibAlign * align = this;
1598   TMatrixD * arrayAlign[72];
1599   TMatrixD * arrayAlignDiff[72];
1600   //
1601   for (Int_t i=0;i<72; i++) {
1602     TMatrixD * mat = new TMatrixD(4,4);
1603     mat->UnitMatrix();
1604     arrayAlign[i]=mat;
1605     arrayAlignDiff[i]=(TMatrixD*)(mat->Clone());
1606   }
1607
1608   TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root");
1609   for (Int_t iter=0; iter<niter;iter++){
1610     printf("Iter=\t%d\n",iter);
1611     for (Int_t is0=0;is0<72; is0++) {
1612       //
1613       //TMatrixD  *mati0 = arrayAlign[is0];
1614       TMatrixD matDiff(4,4);
1615       Double_t sumw=0;
1616       for (Int_t is1=0;is1<72; is1++) {
1617         Bool_t invers=kFALSE;
1618         Int_t npoints=0;
1619         TMatrixD covar;
1620         TVectorD errors;
1621         const TMatrixD *mat = align->GetTransformation(is0,is1,0);
1622         if (mat){
1623           npoints = align->GetFitter6(is0,is1)->GetNpoints();
1624           if (npoints>minPoints){
1625             align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar);
1626             align->GetFitter6(is0,is1)->GetErrors(errors);
1627           }
1628         }
1629         else{
1630           invers=kTRUE;
1631           mat = align->GetTransformation(is1,is0,0);
1632           if (mat) {
1633             npoints = align->GetFitter6(is1,is0)->GetNpoints();
1634             if (npoints>minPoints){
1635               align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar);
1636               align->GetFitter6(is1,is0)->GetErrors(errors);
1637             }
1638           }
1639         }
1640         if (!mat) continue;
1641         if (npoints<minPoints) continue;
1642         //
1643         Double_t weight=1;
1644         if (is1/36>is0/36) weight*=2./3.; //IROC-OROC
1645         if (is1/36<is0/36) weight*=1./3.; //OROC-IROC
1646         if (is1/36==is0/36) weight*=1/3.; //OROC-OROC
1647         if (is1%36!=is0%36) weight*=1/2.; //Not up-down
1648         weight/=(errors[4]*errors[4]+sysError*sysError);
1649         //
1650         //
1651         TMatrixD matT = *mat;
1652         if (invers)  matT.Invert();
1653         TMatrixD diffMat= (*(arrayAlign[is1]))*matT;
1654         diffMat-=(*arrayAlign[is0]);
1655         matDiff+=weight*diffMat;
1656         sumw+=weight;
1657
1658         (*cstream)<<"LAlign"<<
1659           "iter="<<iter<<
1660           "s0="<<is0<<
1661           "s1="<<is1<<
1662           "npoints="<<npoints<<
1663           "m60.="<<arrayAlign[is0]<<
1664           "m61.="<<arrayAlign[is1]<<
1665           "m01.="<<&matT<<
1666           "diff.="<<&diffMat<<
1667           "cov.="<<&covar<<
1668           "err.="<<&errors<<
1669           "\n";
1670       }
1671       if (sumw>0){
1672         matDiff*=1/sumw;
1673         matDiff(0,0)=0;
1674         matDiff(1,1)=0;
1675         matDiff(1,1)=0;
1676         matDiff(1,1)=0;
1677         (*arrayAlignDiff[is0]) = matDiff;
1678       }
1679     }
1680     for (Int_t is0=0;is0<72; is0++) {
1681       if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]);
1682       if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]);
1683        //
1684       (*cstream)<<"GAlign"<<
1685         "iter="<<iter<<
1686         "s0="<<is0<<
1687         "m6.="<<arrayAlign[is0]<<
1688         "\n";
1689     }
1690   }
1691   delete cstream;
1692   for (Int_t isec=0;isec<72;isec++){
1693     fCombinedMatrixArray6.AddAt(arrayAlign[isec],isec);
1694     delete arrayAlignDiff[isec];
1695   }
1696 }
1697
1698
1699 /*
1700
1701
1702 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1703 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1704 AliXRDPROOFtoolkit tool;
1705 TChain * chainTr = tool.MakeChain("align.txt","Track",0,10200);
1706 chainTr->Lookup();
1707 TChain * chainTracklet = tool.MakeChain("align.txt","Tracklet",0,20);
1708 chainTracklet->Lookup();
1709
1710
1711 TCut cutS0("sqrt(tp2.fC[0]+tp2.fC[0])<0.6");
1712 TCut cutS1("sqrt(tp2.fC[2]+tp2.fC[2])<0.6");
1713 TCut cutS2("sqrt(tp2.fC[5]+tp2.fC[5])<0.04");
1714 TCut cutS3("sqrt(tp2.fC[9]+tp2.fC[9])<0.02");
1715 TCut cutS4("sqrt(tp2.fC[14]+tp2.fC[14])<0.5");
1716 // resolution cuts
1717 TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
1718
1719 TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1720 TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1721 TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.02");
1722 TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.01");
1723 TCut cutE("abs(tp2.fP[1])<235");
1724 TCut cutP=cutP0+cutP1+cutP2+cutP3+cutE;
1725
1726
1727 //
1728 //
1729 TCut cutA =  cutP+cutS;
1730 chainTr->Draw(">>listEL",cutA,"entryList");
1731 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1732 chainTr->SetEntryList(elist);
1733
1734
1735
1736
1737
1738
1739 TCut cutS("s1%36==s2%36");
1740
1741 TCut cutN("c1>32&&c2>60");
1742 TCut cutC0("sqrt(tp2.fC[0])<1");
1743
1744 TCut cutX("abs(tp2.fX-133.6)<2");
1745
1746 TCut cutA = cutP+cutN;
1747
1748
1749 TCut cutY("abs(vcY.fElements-vtY.fElements)<0.3&&vcY.fElements!=0")
1750 TCut cutZ("abs(vcZ.fElements-vtZ.fElements)<0.3&&vcZ.fElements!=0")
1751
1752
1753
1754
1755
1756
1757 TCut cutA =  cutP+cutS;
1758 chainTracklet->Draw(">>listEL",cutA,"entryList");
1759 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1760 chainTracklet->SetEntryList(elist);
1761
1762 //
1763 TVectorD * vec1 = 0;
1764 TVectorD * vec2 = 0;
1765 AliExternalTrackParam * tp1 = 0;
1766 AliExternalTrackParam * tp2 = 0;
1767
1768 Int_t      s1 = 0;
1769 Int_t      s2 = 0;
1770 chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1771 chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1772 chainTracklet->GetBranch("s1")->SetAddress(&s1);
1773 chainTracklet->GetBranch("s2")->SetAddress(&s2);
1774
1775
1776 AliTPCcalibAlign align;
1777 {
1778 for (Int_t i=0; i< elist->GetN(); i++){
1779 //for (Int_t i=0; i<100000; i++){
1780 chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
1781 chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
1782 chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1783 chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1784 chainTracklet->GetBranch("s1")->SetAddress(&s1);
1785 chainTracklet->GetBranch("s2")->SetAddress(&s2);
1786
1787 chainTracklet->GetEntry(i);
1788 if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i,tentry, s1,s2);
1789 //vec1.Print();
1790 TLinearFitter * fitter = align.GetOrMakeFitter6(s1,s2);
1791 if (fitter) align.Process6(vec1->GetMatrixArray(),vec2->GetMatrixArray(),fitter);
1792 }
1793 }
1794
1795 */