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