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