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 ////
  72
  73 #include "TLinearFitter.h"
  74 #include "AliTPCcalibAlign.h"
  75 #include "AliExternalTrackParam.h"
  76 #include "AliTPCTracklet.h"
  77 #include "TH1D.h"
  78 #include "TVectorD.h"
  79 #include "TTreeStream.h"
  80 #include "TFile.h"
  81 #include "TF1.h"
  82 #include "TGraphErrors.h"
  83 #include "AliTPCclusterMI.h"
  84 #include "AliTPCseed.h"
  85 #include "AliTracker.h"
  86 #include "TClonesArray.h"
  87
  88
  89 #include "TTreeStream.h"
  90 #include <iostream>
  91 #include <sstream>
  92 using namespace std;
  93
  94 ClassImp(AliTPCcalibAlign)
  95
  96 AliTPCcalibAlign::AliTPCcalibAlign()
  97   :  fDphiHistArray(72*72),
  98      fDthetaHistArray(72*72),
  99      fDyHistArray(72*72),
 100      fDzHistArray(72*72),
 101      fFitterArray12(72*72),
 102      fFitterArray9(72*72),
 103      fFitterArray6(72*72)
 104 {
 105   //
 106   // Constructor
 107   //
 108   for (Int_t i=0;i<72*72;++i) {
 109     fPoints[i]=0;
 110   }
 111 }
 112
 113 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
 114   :AliTPCcalibBase(),
 115    fDphiHistArray(72*72),
 116    fDthetaHistArray(72*72),
 117    fDyHistArray(72*72),
 118    fDzHistArray(72*72),
 119    fFitterArray12(72*72),
 120    fFitterArray9(72*72),
 121    fFitterArray6(72*72)
 122 {
 123   //
 124   // Constructor
 125   //
 126   SetName(name);
 127   SetTitle(title);
 128   for (Int_t i=0;i<72*72;++i) {
 129     fPoints[i]=0;
 130   }
 131 }
 132
 133 AliTPCcalibAlign::~AliTPCcalibAlign() {
 134   //
 135   // destructor
 136   //
 137 }
 138
 139 void AliTPCcalibAlign::Process(AliTPCseed *seed) {
 140   TObjArray tracklets=
 141     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
 142                                     kFALSE,20,2);
 143  //  TObjArray trackletsL=
 144 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
 145 //                                  kFALSE,20,2);
 146 //   TObjArray trackletsQ=
 147 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
 148 //                                  kFALSE,20,2);
 149 //   TObjArray trackletsR=
 150 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
 151 //                                  kFALSE,20,2);
 152   tracklets.SetOwner();
 153  //  trackletsL.SetOwner();
 154 //   trackletsQ.SetOwner();
 155 //   trackletsR.SetOwner();
 156   if (tracklets.GetEntries()==2) {
 157     AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
 158     AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
 159     if (t1->GetSector()>t2->GetSector()) {
 160       AliTPCTracklet* tmp=t1;
 161       t1=t2;
 162       t2=tmp;
 163     }
 164     AliExternalTrackParam *common1=0,*common2=0;
 165     if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
 166       ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
 167     delete common1;
 168     delete common2;
 169   }
 170
 171 }
 172
 173 void AliTPCcalibAlign::Analyze(){
 174   //
 175   // Analyze function
 176   //
 177   EvalFitters();
 178 }
 179
 180
 181 void AliTPCcalibAlign::Terminate(){
 182   //
 183   // Terminate function
 184   // call base terminate + Eval of fitters
 185   //
 186   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Terminate");
 187   EvalFitters();
 188   AliTPCcalibBase::Terminate();
 189 }
 190
 191
 192
 193
 194 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
 195                                         const AliExternalTrackParam &tp2,
 196                                         const AliTPCseed * seed,
 197                                         Int_t s1,Int_t s2) {
 198
 199   //
 200   //
 201   //
 202   //
 203   //
 204   // Process function to fill fitters
 205   //
 206   Double_t t1[5],t2[5];
 207   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 208   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 209   x1   =tp1.GetX();
 210   y1   =tp1.GetY();
 211   z1   =tp1.GetZ();
 212   Double_t snp1=tp1.GetSnp();
 213   dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
 214   Double_t tgl1=tp1.GetTgl();
 215   // dz/dx = 1/(cos(theta)*cos(phi))
 216   dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1);
 217   x2   =tp2.GetX();
 218   y2   =tp2.GetY();
 219   z2   =tp2.GetZ();
 220   Double_t snp2=tp2.GetSnp();
 221   dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
 222   Double_t tgl2=tp2.GetTgl();
 223   dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2);
 224   //
 225   //
 226   //
 227   if (fStreamLevel>1){
 228     TTreeSRedirector *cstream = GetDebugStreamer();
 229     if (cstream){
 230       static TVectorD vec1(5);
 231       static TVectorD vec2(5);
 232       vec1.SetElements(t1);
 233       vec2.SetElements(t2);
 234       AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
 235       AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
 236       (*cstream)<<"Tracklet"<<
 237         "tp1.="<<p1<<
 238         "tp2.="<<p2<<
 239         "v1.="<<&vec1<<
 240         "v2.="<<&vec2<<
 241         "s1="<<s1<<
 242         "s2="<<s2<<
 243         "\n";
 244     }
 245   }
 246   //
 247   // Aplly cut selection
 248   /*
 249     TChain * chainalign = tool.MakeChain("align.txt","Tracklet",0,1000000)
 250     chainalign->Lookup();
 251
 252     // Cuts to be justified with the debug streamer
 253     //
 254     TCut c1pt("abs((tp1.fP[4]+tp2.fP[4])*0.5)<3"); // pt cut  - OK
 255     TCut cdy("abs(tp1.fP[0]-tp2.fP[0])<2");
 256     TCut cdz("abs(tp1.fP[1]-tp2.fP[1])<2");
 257     TCut cdphi("abs(tp1.fP[2]-tp2.fP[2])<0.02");
 258     TCut cdt("abs(tp1.fP[3]-tp2.fP[3])<0.02");
 259     TCut cd1pt("abs(tp1.fP[4]-tp2.fP[4])<0.3");    // delta 1/pt cut  -OK
 260     //
 261     //
 262     TCut acut =  c1pt+cdy+cdz+cdphi+cdt+cd1pt;
 263
 264   */
 265   //   1. pt cut
 266   //   2. dy
 267   //   3. dz
 268   //   4. dphi
 269   //   5. dtheta
 270   //   6. d1pt
 271   printf("Process track\n");
 272   if (TMath::Abs(tp1.GetParameter()[0]-tp2.GetParameter()[0])>2)    return;
 273   if (TMath::Abs(tp1.GetParameter()[1]-tp2.GetParameter()[1])>2)    return;
 274   if (TMath::Abs(tp1.GetParameter()[2]-tp2.GetParameter()[2])>0.02) return;
 275   if (TMath::Abs(tp1.GetParameter()[3]-tp2.GetParameter()[3])>0.02) return;
 276   if (TMath::Abs(tp1.GetParameter()[4]-tp2.GetParameter()[4])>0.3)  return;
 277   if (TMath::Abs((tp1.GetParameter()[4]+tp2.GetParameter()[4])*0.5)>3)  return;
 278   if (TMath::Abs((tp1.GetParameter()[0]-tp2.GetParameter()[0]))<0.000000001)  return;
 279   printf("Filling track\n");
 280   //
 281   // fill resolution histograms - previous cut included
 282   FillHisto(tp1,tp2,s1,s2);
 283   //
 284   Process12(t1,t2,GetOrMakeFitter12(s1,s2));
 285   Process9(t1,t2,GetOrMakeFitter9(s1,s2));
 286   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
 287   ProcessDiff(tp1,tp2, seed,s1,s2);
 288   ++fPoints[GetIndex(s1,s2)];
 289 }
 290
 291 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
 292                                     const AliExternalTrackParam &t2,
 293                                     const AliTPCseed *seed,
 294                                     Int_t s1,Int_t s2)
 295 {
 296   //
 297   // Process local residuals function
 298   //
 299   TVectorD vecX(160);
 300   TVectorD vecY(160);
 301   TVectorD vecZ(160);
 302   TVectorD vecClY(160);
 303   TVectorD vecClZ(160);
 304   TClonesArray arrCl("AliTPCclusterMI",160);
 305   arrCl.ExpandCreateFast(160);
 306   Int_t count1=0, count2=0;
 307   for (Int_t i=0;i<160;++i) {
 308     AliTPCclusterMI *c=seed->GetClusterPointer(i);
 309     vecX[i]=0;
 310     vecY[i]=0;
 311     vecZ[i]=0;
 312     if (!c) continue;
 313     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
 314     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
 315     vecClY[i] = c->GetY();
 316     vecClZ[i] = c->GetZ();
 317     cl=*c;
 318     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
 319     if (c->GetDetector()==s1) ++count1;
 320     if (c->GetDetector()==s2) ++count2;
 321     Double_t gxyz[3],xyz[3];
 322     t1.GetXYZ(gxyz);
 323     Float_t bz = AliTracker::GetBz(gxyz);
 324     par->GetYAt(c->GetX(), bz, xyz[1]);
 325     par->GetZAt(c->GetX(), bz, xyz[2]);
 326     vecX[i] = c->GetX();
 327     vecY[i]= xyz[1];
 328     vecZ[i]= xyz[2];
 329   }
 330   //
 331   //
 332   if (fStreamLevel>5){
 333     //
 334     // huge output - cluster residuals to be investigated
 335     //
 336     TTreeSRedirector *cstream = GetDebugStreamer();
 337     AliTPCseed * t = (AliTPCseed*) seed;
 338     //AliExternalTrackParam *p0 = &((AliExternalTrackParam&)seed);
 339     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
 340     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
 341     /*
 342
 343       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");
 344
 345     */
 346
 347     if (cstream){
 348       (*cstream)<<"Track"<<
 349         "Cl.="<<&arrCl<<
 350         //"tp0.="<<p0<<
 351         "tp1.="<<p1<<
 352         "tp2.="<<p2<<
 353         "vtX.="<<&vecX<<
 354         "vtY.="<<&vecY<<
 355         "vtZ.="<<&vecZ<<
 356         "vcY.="<<&vecClY<<
 357         "vcZ.="<<&vecClZ<<
 358         "s1="<<s1<<
 359         "s2="<<s2<<
 360         "c1="<<count1<<
 361         "c2="<<count2<<
 362         "\n";
 363     }
 364   }
 365 }
 366
 367
 368
 369
 370 void AliTPCcalibAlign::Process12(const Double_t *t1,
 371                                  const Double_t *t2,
 372                                  TLinearFitter *fitter) {
 373   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 374   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 375   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 376   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 377   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 378   //
 379   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
 380   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
 381   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
 382
 383
 384
 385   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 386   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 387
 388   // TODO:
 389   Double_t sy    = 0.1;
 390   Double_t sz    = 0.1;
 391   Double_t sdydx = 0.001;
 392   Double_t sdzdx = 0.001;
 393
 394   Double_t p[12];
 395   Double_t value;
 396
 397   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 398   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 399   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 400   for (Int_t i=0; i<12;i++) p[i]=0.;
 401   p[3+0] = x1;          // a10
 402   p[3+1] = y1;          // a11
 403   p[3+2] = z1;          // a12
 404   p[9+1] = 1.;          // a13
 405   p[0+1] = y1*dydx2;    // a01
 406   p[0+2] = z1*dydx2;    // a02
 407   p[9+0] = dydx2;       // a03
 408   value  = y2;
 409   fitter->AddPoint(p,value,sy);
 410
 411   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 412   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 413   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 414   for (Int_t i=0; i<12;i++) p[i]=0.;
 415   p[6+0] = x1;           // a20
 416   p[6+1] = y1;           // a21
 417   p[6+2] = z1;           // a22
 418   p[9+2] = 1.;           // a23
 419   p[0+1] = y1*dzdx2;     // a01
 420   p[0+2] = z1*dzdx2;     // a02
 421   p[9+0] = dzdx2;        // a03
 422   value  = z2;
 423   fitter->AddPoint(p,value,sz);
 424
 425   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 426   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 427   for (Int_t i=0; i<12;i++) p[i]=0.;
 428   p[3+0] = 1.;           // a10
 429   p[3+1] = dydx1;        // a11
 430   p[3+2] = dzdx1;        // a12
 431   p[0+0] = -dydx2;       // a00
 432   p[0+1] = -dydx1*dydx2; // a01
 433   p[0+2] = -dzdx1*dydx2; // a02
 434   value  = 0.;
 435   fitter->AddPoint(p,value,sdydx);
 436
 437   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 438   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 439   for (Int_t i=0; i<12;i++) p[i]=0.;
 440   p[6+0] = 1;            // a20
 441   p[6+1] = dydx1;        // a21
 442   p[6+2] = dzdx1;        // a22
 443   p[0+0] = -dzdx2;       // a00
 444   p[0+1] = -dydx1*dzdx2; // a01
 445   p[0+2] = -dzdx1*dzdx2; // a02
 446   value  = 0.;
 447   fitter->AddPoint(p,value,sdzdx);
 448 }
 449
 450 void AliTPCcalibAlign::Process9(Double_t *t1,
 451                                 Double_t *t2,
 452                                 TLinearFitter *fitter) {
 453   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 454   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 455   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 456   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 457   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 458   //
 459   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
 460   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
 461   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8]
 462
 463
 464   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 465   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 466
 467   // TODO:
 468   Double_t sy    = 0.1;
 469   Double_t sz    = 0.1;
 470   Double_t sdydx = 0.001;
 471   Double_t sdzdx = 0.001;
 472   //
 473   Double_t p[12];
 474   Double_t value;
 475
 476   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 477   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 478   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 479   for (Int_t i=0; i<12;i++) p[i]=0.;
 480   p[2]   += x1;           // a10
 481   //p[]  +=1;             // a11
 482   p[3]   += z1;           // a12
 483   p[7]   += 1;            // a13
 484   p[0]   += y1*dydx2;     // a01
 485   p[1]   += z1*dydx2;     // a02
 486   p[6]   += dydx2;        // a03
 487   value   = y2-y1;        //-a11
 488   fitter->AddPoint(p,value,sy);
 489   //
 490   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 491   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 492   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 493   for (Int_t i=0; i<12;i++) p[i]=0.;
 494   p[4]   += x1;           // a20
 495   p[5]   += y1;           // a21
 496   //p[]  += z1;           // a22
 497   p[8]   += 1.;           // a23
 498   p[0]   += y1*dzdx2;     // a01
 499   p[1]   += z1*dzdx2;     // a02
 500   p[6]   += dzdx2;        // a03
 501   value  = z2-z1;         //-a22
 502   fitter->AddPoint(p,value,sz);
 503
 504   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 505   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 506   for (Int_t i=0; i<12;i++) p[i]=0.;
 507   p[2]   += 1.;           // a10
 508   //p[]  += dydx1;      // a11
 509   p[3]   += dzdx1;        // a12
 510   //p[]  += -dydx2;       // a00
 511   p[0]   += -dydx1*dydx2; // a01
 512   p[1]   += -dzdx1*dydx2; // a02
 513   value  = -dydx1+dydx2;  // -a11 + a00
 514   fitter->AddPoint(p,value,sdydx);
 515
 516   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 517   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 518   for (Int_t i=0; i<12;i++) p[i]=0.;
 519   p[4]   += 1;            // a20
 520   p[5]   += dydx1;        // a21
 521   //p[]  += dzdx1;        // a22
 522   //p[]  += -dzdx2;       // a00
 523   p[0]   += -dydx1*dzdx2; // a01
 524   p[1]   += -dzdx1*dzdx2; // a02
 525   value  = -dzdx1+dzdx2;  // -a22 + a00
 526   fitter->AddPoint(p,value,sdzdx);
 527 }
 528
 529 void AliTPCcalibAlign::Process6(Double_t *t1,
 530                                 Double_t *t2,
 531                                 TLinearFitter *fitter) {
 532   // x2    =  1  *x1 +-a01*y1 + 0      +a03
 533   // y2    =  a01*x1 + 1  *y1 + 0      +a13
 534   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
 535   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 536   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 537   //
 538   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
 539   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
 540   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5]
 541
 542   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 543   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 544
 545   // TODO:
 546   Double_t sy    = 0.1;
 547   Double_t sz    = 0.1;
 548   Double_t sdydx = 0.001;
 549   Double_t sdzdx = 0.001;
 550
 551   Double_t p[12];
 552   Double_t value;
 553   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 554   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 555   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 556   for (Int_t i=0; i<12;i++) p[i]=0.;
 557   p[0]   += x1;           // a10
 558   //p[]  +=1;             // a11
 559   //p[]  += z1;           // a12
 560   p[4]   += 1;            // a13
 561   p[0]   += -y1*dydx2;    // a01
 562   //p[]  += z1*dydx2;     // a02
 563   p[3]   += dydx2;        // a03
 564   value   = y2-y1;        //-a11
 565   fitter->AddPoint(p,value,sy);
 566   //
 567   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 568   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 569   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 570   for (Int_t i=0; i<12;i++) p[i]=0.;
 571   p[1]   += x1;           // a20
 572   p[2]   += y1;           // a21
 573   //p[]  += z1;           // a22
 574   p[5]   += 1.;           // a23
 575   p[0]   += -y1*dzdx2;    // a01
 576   //p[]   += z1*dzdx2;     // a02
 577   p[3]   += dzdx2;        // a03
 578   value  = z2-z1;         //-a22
 579   fitter->AddPoint(p,value,sz);
 580
 581   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 582   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 583   for (Int_t i=0; i<12;i++) p[i]=0.;
 584   p[0]   += 1.;           // a10
 585   //p[]  += dydx1;      // a11
 586   //p[]   += dzdx1;        // a12
 587   //p[]  += -dydx2;       // a00
 588   p[0]   +=  dydx1*dydx2; // a01
 589   //p[]   += -dzdx1*dydx2; // a02
 590   value  = -dydx1+dydx2;  // -a11 + a00
 591   fitter->AddPoint(p,value,sdydx);
 592
 593   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 594   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 595   for (Int_t i=0; i<12;i++) p[i]=0.;
 596   p[1]   += 1;            // a20
 597   p[2]   += dydx1;        // a21
 598   //p[]  += dzdx1;        // a22
 599   //p[]  += -dzdx2;       // a00
 600   p[0]   +=  dydx1*dzdx2; // a01
 601   //p[]  += -dzdx1*dzdx2; // a02
 602   value  = -dzdx1+dzdx2;  // -a22 + a00
 603   fitter->AddPoint(p,value,sdzdx);
 604 }
 605
 606
 607
 608
 609 void AliTPCcalibAlign::EvalFitters() {
 610   //
 611   // Analyze function
 612   //
 613   // Perform the fitting using linear fitters
 614   //
 615   Int_t kMinPoints =50;
 616   TLinearFitter *f;
 617   TFile fff("alignDebug.root","recreate");
 618   for (Int_t s1=0;s1<72;++s1)
 619     for (Int_t s2=0;s2<72;++s2){
 620       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 621         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 622         if (f->Eval()!=0) {
 623           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 624           f->Write(Form("f12_%d_%d",s1,s2));
 625         }else{
 626           f->Write(Form("f12_%d_%d",s1,s2));
 627         }
 628       }
 629       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 630         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 631         if (f->Eval()!=0) {
 632           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 633         }else{
 634           f->Write(Form("f9_%d_%d",s1,s2));
 635         }
 636       }
 637       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 638         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 639         if (f->Eval()!=0) {
 640           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 641         }else{
 642           f->Write(Form("f6_%d_%d",s1,s2));
 643         }
 644       }
 645     }
 646   this->Write("align");
 647   /*
 648
 649   fitter->Eval();
 650   fitter->Eval();
 651   chi212 = align->GetChisquare()/(4.*entries);
 652
 653   TMatrixD mat(13,13);
 654   TVectorD par(13);
 655   align->GetParameters(par);
 656   align->GetCovarianceMatrix(mat);
 657
 658   //
 659   //
 660   for (Int_t i=0; i<12;i++){
 661     palign12(i)= par(i+1);
 662     for (Int_t j=0; j<12;j++){
 663       pcovar12(i,j)   = mat(i+1,j+1);
 664       pcovar12(i,j) *= chi212;
 665     }
 666   }
 667   //
 668   for (Int_t i=0; i<12;i++){
 669     psigma12(i)  = TMath::Sqrt(pcovar12(i,i));
 670     palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i));
 671     for (Int_t j=0; j<12;j++){
 672       pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j));
 673     }
 674   }
 675   */
 676 }
 677
 678 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
 679   //
 680   // get or make fitter - general linear transformation
 681   //
 682   static Int_t counter12=0;
 683   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]");
 684   TLinearFitter * fitter = GetFitter12(s1,s2);
 685   if (fitter) return fitter;
 686   //  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]");
 687   fitter =new TLinearFitter(&f12,"");
 688   fitter->StoreData(kFALSE);
 689   fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
 690   counter12++;
 691   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
 692   return fitter;
 693 }
 694
 695 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
 696   //
 697   //get or make fitter - general linear transformation - no scaling
 698   //
 699   static Int_t counter9=0;
 700   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 701   TLinearFitter * fitter = GetFitter9(s1,s2);
 702   if (fitter) return fitter;
 703   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 704   fitter =new TLinearFitter(&f9,"");
 705   fitter->StoreData(kFALSE);
 706   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
 707   counter9++;
 708   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
 709   return fitter;
 710 }
 711
 712 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
 713   //
 714   // get or make fitter  - 6 paramater linear tranformation
 715   //                     - no scaling
 716   //                     - rotation x-y
 717   //                     - tilting x-z, y-z
 718   static Int_t counter6=0;
 719   static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
 720   TLinearFitter * fitter = GetFitter6(s1,s2);
 721   if (fitter) return fitter;
 722   //  fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
 723   fitter=new TLinearFitter(&f6,"");
 724   fitter->StoreData(kFALSE);
 725   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
 726   counter6++;
 727   if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
 728   return fitter;
 729 }
 730
 731
 732
 733
 734
 735 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
 736   //
 737   // GetTransformation matrix - 12 paramaters - generael linear transformation
 738   //
 739   if (!GetFitter12(s1,s2))
 740     return false;
 741   else {
 742     TVectorD p(12);
 743     GetFitter12(s1,s2)->GetParameters(p);
 744     a.ResizeTo(4,4);
 745     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
 746     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
 747     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
 748     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 749     return true;
 750   }
 751 }
 752
 753 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
 754   //
 755   // GetTransformation matrix - 9 paramaters - general linear transformation
 756   //                            No scaling
 757   //
 758   if (!GetFitter9(s1,s2))
 759     return false;
 760   else {
 761     TVectorD p(9);
 762     GetFitter9(s1,s2)->GetParameters(p);
 763     a.ResizeTo(4,4);
 764     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
 765     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
 766     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
 767     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 768     return true;
 769   }
 770 }
 771
 772 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
 773   //
 774   // GetTransformation matrix - 6  paramaters
 775   //                            3  translation
 776   //                            1  rotation -x-y
 777   //                            2  tilting x-z y-z
 778   if (!GetFitter6(s1,s2))
 779     return false;
 780   else {
 781     TVectorD p(6);
 782     GetFitter6(s1,s2)->GetParameters(p);
 783     a.ResizeTo(4,4);
 784     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
 785     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
 786     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
 787     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 788     return true;
 789   }
 790 }
 791
 792 void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
 793                                         const AliExternalTrackParam &tp2,
 794                                         Int_t s1,Int_t s2) {
 795   //
 796   // Fill residual histograms
 797   // Innner-Outer
 798   // Left right - x-y
 799   // A-C side
 800   if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0)  {
 801     GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
 802     GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
 803     GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
 804     GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
 805   }
 806 }
 807
 808
 809
 810 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
 811 {
 812   //
 813   // return specified residual histogram - it is only QA
 814   // if force specified the histogram and given histogram is not existing
 815   //  new histogram is created
 816   //
 817   if (GetIndex(s1,s2)>=72*72) return 0;
 818   TObjArray *histoArray=0;
 819   switch (type) {
 820   case kY:
 821     histoArray = &fDyHistArray; break;
 822   case kZ:
 823     histoArray = &fDzHistArray; break;
 824   case kPhi:
 825     histoArray = &fDphiHistArray; break;
 826   case kTheta:
 827     histoArray = &fDthetaHistArray; break;
 828   }
 829   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
 830   if (histo) return histo;
 831   if (force==kFALSE) return 0;
 832   //
 833   stringstream name;
 834   stringstream title;
 835   switch (type) {
 836   case kY:
 837     name<<"hist_y_"<<s1<<"_"<<s2;
 838     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
 839     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
 840     break;
 841   case kZ:
 842     name<<"hist_z_"<<s1<<"_"<<s2;
 843     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
 844     histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
 845     break;
 846   case kPhi:
 847     name<<"hist_phi_"<<s1<<"_"<<s2;
 848     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
 849     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
 850     break;
 851   case kTheta:
 852     name<<"hist_theta_"<<s1<<"_"<<s2;
 853     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
 854     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
 855     break;
 856   }
 857   histo->SetDirectory(0);
 858   histoArray->AddAt(histo,GetIndex(s1,s2));
 859   return histo;
 860 }
 861
 862 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
 863                                            Int_t i0, Int_t i1, FitType type)
 864 {
 865   //
 866   //
 867   //
 868   TMatrixD mat;
 869   TObjArray *fitArray=0;
 870   Double_t xsec[1000];
 871   Double_t ysec[1000];
 872   Int_t npoints=0;
 873   for (Int_t isec = sec0; isec<=sec1; isec++){
 874     Int_t isec2 = (isec+dsec)%72;
 875     switch (type) {
 876     case k6:
 877       GetTransformation6(isec,isec2,mat);break;
 878     case k9:
 879       GetTransformation9(isec,isec2,mat);break;
 880     case k12:
 881       GetTransformation12(isec,isec2,mat);break;
 882     }
 883     xsec[npoints]=isec;
 884     ysec[npoints]=mat(i0,i1);
 885     ++npoints;
 886   }
 887   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
 888   Char_t name[1000];
 889   sprintf(name,"Mat[%d,%d]  Type=%d",i0,i1,type);
 890   gr->SetName(name);
 891   return gr;
 892 }
 893
 894 void  AliTPCcalibAlign::MakeTree(const char *fname){
 895   //
 896   // make tree with alignment cosntant  -
 897   // For  QA visualization
 898   //
 899   /*
 900    TFile f("CalibObjects.root");
 901    TObjArray *array  = (TObjArray*)f.Get("TPCCalib");
 902    AliTPCcalibAlign   *alignTPC = (AliTPCcalibAlign   *)array->At(0);
 903    alignTPC->MakeTree("alignTree.root");
 904    TFile falign("alignTree.root");
 905    Align->Draw("dy")
 906    */
 907   const Int_t kMinPoints=50;
 908   TTreeSRedirector cstream(fname);
 909   for (Int_t s1=0;s1<72;++s1)
 910     for (Int_t s2=0;s2<72;++s2){
 911       if (fPoints[GetIndex(s1,s2)]<kMinPoints) continue;
 912       TMatrixD m6;
 913       TMatrixD m9;
 914       TMatrixD m12;
 915       GetTransformation6(s1,s2,m6);
 916       GetTransformation9(s1,s2,m9);
 917       GetTransformation12(s1,s2,m12);
 918       Double_t dy=0, dz=0, dphi=0,dtheta=0;
 919       Double_t sy=0, sz=0, sphi=0,stheta=0;
 920       Double_t ny=0, nz=0, nphi=0,ntheta=0;
 921       TH1 * his=0;
 922       his = GetHisto(kY,s1,s2);
 923       if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
 924       his = GetHisto(kZ,s1,s2);
 925       if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
 926       his = GetHisto(kPhi,s1,s2);
 927       if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
 928       his = GetHisto(kTheta,s1,s2);
 929       if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
 930       //
 931       cstream<<"Align"<<
 932         "s1="<<s1<<     // reference sector
 933         "s2="<<s2<<     // sector to align
 934         "m6.="<<&m6<<   // tranformation matrix
 935         "m9.="<<&m9<<   //
 936         "m12.="<<&m12<<
 937         //               histograms mean RMS and entries
 938         "dy="<<dy<<
 939         "sy="<<sy<<
 940         "ny="<<ny<<
 941         "dz="<<dz<<
 942         "sz="<<sz<<
 943         "nz="<<nz<<
 944         "dphi="<<dphi<<
 945         "sphi="<<sphi<<
 946         "nphi="<<nphi<<
 947         "dtheta="<<dtheta<<
 948         "stheta="<<stheta<<
 949         "ntheta="<<ntheta<<
 950         "\n";
 951     }
 952
 953 }
 954
 955
 956 //_____________________________________________________________________
 957 Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
 958   //
 959   // merge function
 960   //
 961   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
 962   if (!list)
 963     return 0;
 964   if (list->IsEmpty())
 965     return 1;
 966
 967   TIterator* iter = list->MakeIterator();
 968   TObject* obj = 0;
 969   iter->Reset();
 970   Int_t count=0;
 971   while((obj = iter->Next()) != 0)
 972     {
 973       AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
 974       if (entry == 0) continue;
 975       Add(entry);
 976       count++;
 977     }
 978   return count;
 979 }
 980
 981
 982 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
 983   //
 984   // Add entry
 985   //
 986
 987   for (Int_t i=0; i<72;i++){
 988     for (Int_t j=0; j<72;j++){
 989       fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
 990
 991       //
 992       // dy
 993       TH1* hdy0 = GetHisto(kY,i,j);
 994       TH1* hdy1 = align->GetHisto(kY,i,j);
 995       if (hdy1){
 996         if (hdy0) hdy0->Add(hdy1);
 997         else {
 998           hdy0 = GetHisto(kY,i,j,kTRUE);
 999           hdy0->Add(hdy1);
1000         }
1001       }
1002       //
1003       // dz
1004       TH1* hdz0 = GetHisto(kZ,i,j);
1005       TH1* hdz1 = align->GetHisto(kZ,i,j);
1006       if (hdz1){
1007         if (hdz0) hdz0->Add(hdz1);
1008         else {
1009           hdz0 = GetHisto(kZ,i,j,kTRUE);
1010           hdz0->Add(hdz1);
1011         }
1012       }
1013       //
1014       // dphi
1015       TH1* hdphi0 = GetHisto(kPhi,i,j);
1016       TH1* hdphi1 = align->GetHisto(kPhi,i,j);
1017       if (hdphi1){
1018         if (hdphi0) hdphi0->Add(hdphi1);
1019         else {
1020           hdphi0 = GetHisto(kPhi,i,j,kTRUE);
1021           hdphi0->Add(hdphi1);
1022         }
1023       }
1024       //
1025       // dtheta
1026       TH1* hdTheta0 = GetHisto(kTheta,i,j);
1027       TH1* hdTheta1 = align->GetHisto(kTheta,i,j);
1028       if (hdTheta1){
1029         if (hdTheta0) hdTheta0->Add(hdTheta1);
1030         else {
1031           hdTheta0 = GetHisto(kTheta,i,j,kTRUE);
1032           hdTheta0->Add(hdTheta1);
1033         }
1034       }
1035     }
1036   }
1037   TLinearFitter *f0=0;
1038   TLinearFitter *f1=0;
1039   for (Int_t i=0; i<72;i++){
1040     for (Int_t j=0; j<72;j++){
1041       //
1042       // fitter12
1043       f0 =  GetFitter12(i,j);
1044       f1 =  GetFitter12(i,j);
1045       if (f1){
1046         if (f0) f0->Add(f1);
1047         else {
1048           f0 = GetOrMakeFitter12(i,j);
1049           f0->Add(f1);
1050         }
1051       }
1052       //
1053       // fitter9
1054       f0 =  GetFitter9(i,j);
1055       f1 =  GetFitter9(i,j);
1056       if (f1){
1057         if (f0) f0->Add(f1);
1058         else {
1059           f0 = GetOrMakeFitter9(i,j);
1060           f0->Add(f1);
1061         }
1062       }
1063       f0 =  GetFitter6(i,j);
1064       f1 =  GetFitter6(i,j);
1065       if (f1){
1066         if (f0) f0->Add(f1);
1067         else {
1068           f0 = GetOrMakeFitter6(i,j);
1069           f0->Add(f1);
1070         }
1071       }
1072     }
1073   }
1074 }