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     // Cuts to be justified with the debug streamer
 250     //
 251     TCut c1pt("abs((tp1.fP[4]+tp2.fP[4])*0.5)<5"); // pt cut
 252     TCut cd1pt("abs(tp1.fP[4]-tp2.fP[4])<0.5");
 253     TCut c1ptpull("abs((tp1.fP[4]-tp2.fP[4])/sqrt(tp1.fC[14]+tp2.fC[14]))<5");
 254     TCut csy("sqrt(tp1.fC[0]+tp2.fC[0])<0.5");
 255     TCut csz("sqrt(tp1.fC[2]+tp2.fC[2])<0.3");
 256     TCut cphi("abs(v1.fElements[3])<1")
 257     TCut ctheta("abs(v1.fElements[4])<1")
 258     //
 259     //
 260     TCut acut =  c1ptpull+c1pt+cd1pt+csy+csz+cphi+ctheta;
 261   */
 262   //   1. pt cut
 263   //   2. delta in curvature
 264   //   3. pull in 1pt
 265   //   4. sigma y
 266   //   5. sigma z
 267   //   6. angle phi
 268   //   7. angle theta
 269   Double_t sigma1pt  = TMath::Sqrt(tp1.GetSigma1Pt2()+tp1.GetSigma1Pt2());
 270   Double_t delta1pt = (tp1.GetParameter()[4]-tp2.GetParameter()[4]);
 271   Double_t pull1pt  = delta1pt/sigma1pt;
 272   if (0.5*TMath::Abs(tp1.GetParameter()[4]+tp2.GetParameter()[4])>5) return;
 273   if (TMath::Abs(delta1pt)>0.5) return;
 274   if (TMath::Abs(pull1pt)>5)    return;
 275   if (TMath::Sqrt(tp1.GetSigmaY2()+tp2.GetSigmaY2())>0.5)    return;
 276   if (TMath::Sqrt(tp1.GetSigmaZ2()+tp2.GetSigmaZ2())>0.3)    return;
 277   if (TMath::Abs(dydx1)>1.)    return;
 278   if (TMath::Abs(dzdx1)>1.)    return;
 279   //
 280   // fill resolution histograms - previous cut included
 281   FillHisto(tp1,tp2,s1,s2);
 282   //
 283   Process12(t1,t2,GetOrMakeFitter12(s1,s2));
 284   Process9(t1,t2,GetOrMakeFitter9(s1,s2));
 285   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
 286   ProcessDiff(tp1,tp2, seed,s1,s2);
 287   ++fPoints[GetIndex(s1,s2)];
 288 }
 289
 290 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
 291                                     const AliExternalTrackParam &t2,
 292                                     const AliTPCseed *seed,
 293                                     Int_t s1,Int_t s2)
 294 {
 295   //
 296   // Process local residuals function
 297   //
 298   TVectorD vecX(160);
 299   TVectorD vecY(160);
 300   TVectorD vecZ(160);
 301   TVectorD vecClY(160);
 302   TVectorD vecClZ(160);
 303   TClonesArray arrCl("AliTPCclusterMI",160);
 304   arrCl.ExpandCreateFast(160);
 305   Int_t count1=0, count2=0;
 306   for (Int_t i=0;i<160;++i) {
 307     AliTPCclusterMI *c=seed->GetClusterPointer(i);
 308     vecX[i]=0;
 309     vecY[i]=0;
 310     vecZ[i]=0;
 311     if (!c) continue;
 312     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
 313     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
 314     vecClY[i] = c->GetY();
 315     vecClZ[i] = c->GetZ();
 316     cl=*c;
 317     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
 318     if (c->GetDetector()==s1) ++count1;
 319     if (c->GetDetector()==s2) ++count2;
 320     Double_t gxyz[3],xyz[3];
 321     t1.GetXYZ(gxyz);
 322     Float_t bz = AliTracker::GetBz(gxyz);
 323     par->GetYAt(c->GetX(), bz, xyz[1]);
 324     par->GetZAt(c->GetX(), bz, xyz[2]);
 325     vecX[i] = c->GetX();
 326     vecY[i]= xyz[1];
 327     vecZ[i]= xyz[2];
 328   }
 329   //
 330   //
 331   if (fStreamLevel>5){
 332     //
 333     // huge output - cluster residuals to be investigated
 334     //
 335     TTreeSRedirector *cstream = GetDebugStreamer();
 336     AliTPCseed * t = (AliTPCseed*) seed;
 337     //AliExternalTrackParam *p0 = &((AliExternalTrackParam&)seed);
 338     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
 339     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
 340     /*
 341
 342       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");
 343
 344     */
 345
 346     if (cstream){
 347       (*cstream)<<"Track"<<
 348         "Cl.="<<&arrCl<<
 349         //"tp0.="<<p0<<
 350         "tp1.="<<p1<<
 351         "tp2.="<<p2<<
 352         "vtX.="<<&vecX<<
 353         "vtY.="<<&vecY<<
 354         "vtZ.="<<&vecZ<<
 355         "vcY.="<<&vecClY<<
 356         "vcZ.="<<&vecClZ<<
 357         "s1="<<s1<<
 358         "s2="<<s2<<
 359         "c1="<<count1<<
 360         "c2="<<count2<<
 361         "\n";
 362     }
 363   }
 364 }
 365
 366
 367
 368
 369 void AliTPCcalibAlign::Process12(const Double_t *t1,
 370                                  const Double_t *t2,
 371                                  TLinearFitter *fitter) {
 372   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 373   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 374   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 375   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 376   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 377   //
 378   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
 379   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
 380   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
 381
 382
 383
 384   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 385   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 386
 387   // TODO:
 388   Double_t sy    = 0.1;
 389   Double_t sz    = 0.1;
 390   Double_t sdydx = 0.001;
 391   Double_t sdzdx = 0.001;
 392
 393   Double_t p[12];
 394   Double_t value;
 395
 396   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 397   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 398   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 399   for (Int_t i=0; i<12;i++) p[i]=0.;
 400   p[3+0] = x1;          // a10
 401   p[3+1] = y1;          // a11
 402   p[3+2] = z1;          // a12
 403   p[9+1] = 1.;          // a13
 404   p[0+1] = y1*dydx2;    // a01
 405   p[0+2] = z1*dydx2;    // a02
 406   p[9+0] = dydx2;       // a03
 407   value  = y2;
 408   fitter->AddPoint(p,value,sy);
 409
 410   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 411   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 412   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 413   for (Int_t i=0; i<12;i++) p[i]=0.;
 414   p[6+0] = x1;           // a20
 415   p[6+1] = y1;           // a21
 416   p[6+2] = z1;           // a22
 417   p[9+2] = 1.;           // a23
 418   p[0+1] = y1*dzdx2;     // a01
 419   p[0+2] = z1*dzdx2;     // a02
 420   p[9+0] = dzdx2;        // a03
 421   value  = z2;
 422   fitter->AddPoint(p,value,sz);
 423
 424   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 425   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 426   for (Int_t i=0; i<12;i++) p[i]=0.;
 427   p[3+0] = 1.;           // a10
 428   p[3+1] = dydx1;        // a11
 429   p[3+2] = dzdx1;        // a12
 430   p[0+0] = -dydx2;       // a00
 431   p[0+1] = -dydx1*dydx2; // a01
 432   p[0+2] = -dzdx1*dydx2; // a02
 433   value  = 0.;
 434   fitter->AddPoint(p,value,sdydx);
 435
 436   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 437   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 438   for (Int_t i=0; i<12;i++) p[i]=0.;
 439   p[6+0] = 1;            // a20
 440   p[6+1] = dydx1;        // a21
 441   p[6+2] = dzdx1;        // a22
 442   p[0+0] = -dzdx2;       // a00
 443   p[0+1] = -dydx1*dzdx2; // a01
 444   p[0+2] = -dzdx1*dzdx2; // a02
 445   value  = 0.;
 446   fitter->AddPoint(p,value,sdzdx);
 447 }
 448
 449 void AliTPCcalibAlign::Process9(Double_t *t1,
 450                                 Double_t *t2,
 451                                 TLinearFitter *fitter) {
 452   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 453   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 454   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 455   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 456   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 457   //
 458   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
 459   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
 460   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8]
 461
 462
 463   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 464   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 465
 466   // TODO:
 467   Double_t sy    = 0.1;
 468   Double_t sz    = 0.1;
 469   Double_t sdydx = 0.001;
 470   Double_t sdzdx = 0.001;
 471   //
 472   Double_t p[12];
 473   Double_t value;
 474
 475   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 476   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 477   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 478   for (Int_t i=0; i<12;i++) p[i]=0.;
 479   p[2]   += x1;           // a10
 480   //p[]  +=1;             // a11
 481   p[3]   += z1;           // a12
 482   p[7]   += 1;            // a13
 483   p[0]   += y1*dydx2;     // a01
 484   p[1]   += z1*dydx2;     // a02
 485   p[6]   += dydx2;        // a03
 486   value   = y2-y1;        //-a11
 487   fitter->AddPoint(p,value,sy);
 488   //
 489   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 490   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 491   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 492   for (Int_t i=0; i<12;i++) p[i]=0.;
 493   p[4]   += x1;           // a20
 494   p[5]   += y1;           // a21
 495   //p[]  += z1;           // a22
 496   p[8]   += 1.;           // a23
 497   p[0]   += y1*dzdx2;     // a01
 498   p[1]   += z1*dzdx2;     // a02
 499   p[6]   += dzdx2;        // a03
 500   value  = z2-z1;         //-a22
 501   fitter->AddPoint(p,value,sz);
 502
 503   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 504   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 505   for (Int_t i=0; i<12;i++) p[i]=0.;
 506   p[2]   += 1.;           // a10
 507   //p[]  += dydx1;      // a11
 508   p[3]   += dzdx1;        // a12
 509   //p[]  += -dydx2;       // a00
 510   p[0]   += -dydx1*dydx2; // a01
 511   p[1]   += -dzdx1*dydx2; // a02
 512   value  = -dydx1+dydx2;  // -a11 + a00
 513   fitter->AddPoint(p,value,sdydx);
 514
 515   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 516   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 517   for (Int_t i=0; i<12;i++) p[i]=0.;
 518   p[4]   += 1;            // a20
 519   p[5]   += dydx1;        // a21
 520   //p[]  += dzdx1;        // a22
 521   //p[]  += -dzdx2;       // a00
 522   p[0]   += -dydx1*dzdx2; // a01
 523   p[1]   += -dzdx1*dzdx2; // a02
 524   value  = -dzdx1+dzdx2;  // -a22 + a00
 525   fitter->AddPoint(p,value,sdzdx);
 526 }
 527
 528 void AliTPCcalibAlign::Process6(Double_t *t1,
 529                                 Double_t *t2,
 530                                 TLinearFitter *fitter) {
 531   // x2    =  1  *x1 +-a01*y1 + 0      +a03
 532   // y2    =  a01*x1 + 1  *y1 + 0      +a13
 533   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
 534   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 535   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 536   //
 537   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
 538   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
 539   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5]
 540
 541   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 542   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 543
 544   // TODO:
 545   Double_t sy    = 0.1;
 546   Double_t sz    = 0.1;
 547   Double_t sdydx = 0.001;
 548   Double_t sdzdx = 0.001;
 549
 550   Double_t p[12];
 551   Double_t value;
 552   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 553   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 554   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 555   for (Int_t i=0; i<12;i++) p[i]=0.;
 556   p[0]   += x1;           // a10
 557   //p[]  +=1;             // a11
 558   //p[]  += z1;           // a12
 559   p[4]   += 1;            // a13
 560   p[0]   += -y1*dydx2;    // a01
 561   //p[]  += z1*dydx2;     // a02
 562   p[3]   += dydx2;        // a03
 563   value   = y2-y1;        //-a11
 564   fitter->AddPoint(p,value,sy);
 565   //
 566   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 567   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 568   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 569   for (Int_t i=0; i<12;i++) p[i]=0.;
 570   p[1]   += x1;           // a20
 571   p[2]   += y1;           // a21
 572   //p[]  += z1;           // a22
 573   p[5]   += 1.;           // a23
 574   p[0]   += -y1*dzdx2;    // a01
 575   //p[]   += z1*dzdx2;     // a02
 576   p[3]   += dzdx2;        // a03
 577   value  = z2-z1;         //-a22
 578   fitter->AddPoint(p,value,sz);
 579
 580   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 581   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 582   for (Int_t i=0; i<12;i++) p[i]=0.;
 583   p[0]   += 1.;           // a10
 584   //p[]  += dydx1;      // a11
 585   //p[]   += dzdx1;        // a12
 586   //p[]  += -dydx2;       // a00
 587   p[0]   +=  dydx1*dydx2; // a01
 588   //p[]   += -dzdx1*dydx2; // a02
 589   value  = -dydx1+dydx2;  // -a11 + a00
 590   fitter->AddPoint(p,value,sdydx);
 591
 592   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 593   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 594   for (Int_t i=0; i<12;i++) p[i]=0.;
 595   p[1]   += 1;            // a20
 596   p[2]   += dydx1;        // a21
 597   //p[]  += dzdx1;        // a22
 598   //p[]  += -dzdx2;       // a00
 599   p[0]   +=  dydx1*dzdx2; // a01
 600   //p[]  += -dzdx1*dzdx2; // a02
 601   value  = -dzdx1+dzdx2;  // -a22 + a00
 602   fitter->AddPoint(p,value,sdzdx);
 603 }
 604
 605
 606
 607
 608 void AliTPCcalibAlign::EvalFitters() {
 609   //
 610   // Analyze function
 611   //
 612   // Perform the fitting using linear fitters
 613   //
 614   Int_t kMinPoints =50;
 615   TLinearFitter *f;
 616   TFile fff("alignDebug.root","recreate");
 617   for (Int_t s1=0;s1<72;++s1)
 618     for (Int_t s2=0;s2<72;++s2){
 619       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 620         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 621         if (f->Eval()!=0) {
 622           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 623           f->Write(Form("f12_%d_%d",s1,s2));
 624         }else{
 625           f->Write(Form("f12_%d_%d",s1,s2));
 626         }
 627       }
 628       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 629         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 630         if (f->Eval()!=0) {
 631           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 632         }else{
 633           f->Write(Form("f9_%d_%d",s1,s2));
 634         }
 635       }
 636       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 637         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 638         if (f->Eval()!=0) {
 639           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 640         }else{
 641           f->Write(Form("f6_%d_%d",s1,s2));
 642         }
 643       }
 644     }
 645   this->Write("align");
 646   /*
 647
 648   fitter->Eval();
 649   fitter->Eval();
 650   chi212 = align->GetChisquare()/(4.*entries);
 651
 652   TMatrixD mat(13,13);
 653   TVectorD par(13);
 654   align->GetParameters(par);
 655   align->GetCovarianceMatrix(mat);
 656
 657   //
 658   //
 659   for (Int_t i=0; i<12;i++){
 660     palign12(i)= par(i+1);
 661     for (Int_t j=0; j<12;j++){
 662       pcovar12(i,j)   = mat(i+1,j+1);
 663       pcovar12(i,j) *= chi212;
 664     }
 665   }
 666   //
 667   for (Int_t i=0; i<12;i++){
 668     psigma12(i)  = TMath::Sqrt(pcovar12(i,i));
 669     palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i));
 670     for (Int_t j=0; j<12;j++){
 671       pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j));
 672     }
 673   }
 674   */
 675 }
 676
 677 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
 678   //
 679   // get or make fitter - general linear transformation
 680   //
 681   static Int_t counter12=0;
 682   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]");
 683   TLinearFitter * fitter = GetFitter12(s1,s2);
 684   if (fitter) return fitter;
 685   //  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]");
 686   fitter =new TLinearFitter(&f12,"");
 687   fitter->StoreData(kFALSE);
 688   fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
 689   counter12++;
 690   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
 691   return fitter;
 692 }
 693
 694 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
 695   //
 696   //get or make fitter - general linear transformation - no scaling
 697   //
 698   static Int_t counter9=0;
 699   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 700   TLinearFitter * fitter = GetFitter9(s1,s2);
 701   if (fitter) return fitter;
 702   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 703   fitter =new TLinearFitter(&f9,"");
 704   fitter->StoreData(kFALSE);
 705   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
 706   counter9++;
 707   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
 708   return fitter;
 709 }
 710
 711 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
 712   //
 713   // get or make fitter  - 6 paramater linear tranformation
 714   //                     - no scaling
 715   //                     - rotation x-y
 716   //                     - tilting x-z, y-z
 717   static Int_t counter6=0;
 718   static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
 719   TLinearFitter * fitter = GetFitter6(s1,s2);
 720   if (fitter) return fitter;
 721   //  fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
 722   fitter=new TLinearFitter(&f6,"");
 723   fitter->StoreData(kFALSE);
 724   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
 725   counter6++;
 726   if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
 727   return fitter;
 728 }
 729
 730
 731
 732
 733
 734 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
 735   //
 736   // GetTransformation matrix - 12 paramaters - generael linear transformation
 737   //
 738   if (!GetFitter12(s1,s2))
 739     return false;
 740   else {
 741     TVectorD p(12);
 742     GetFitter12(s1,s2)->GetParameters(p);
 743     a.ResizeTo(4,4);
 744     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
 745     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
 746     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
 747     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 748     return true;
 749   }
 750 }
 751
 752 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
 753   //
 754   // GetTransformation matrix - 9 paramaters - general linear transformation
 755   //                            No scaling
 756   //
 757   if (!GetFitter9(s1,s2))
 758     return false;
 759   else {
 760     TVectorD p(9);
 761     GetFitter9(s1,s2)->GetParameters(p);
 762     a.ResizeTo(4,4);
 763     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
 764     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
 765     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
 766     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 767     return true;
 768   }
 769 }
 770
 771 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
 772   //
 773   // GetTransformation matrix - 6  paramaters
 774   //                            3  translation
 775   //                            1  rotation -x-y
 776   //                            2  tilting x-z y-z
 777   if (!GetFitter6(s1,s2))
 778     return false;
 779   else {
 780     TVectorD p(6);
 781     GetFitter6(s1,s2)->GetParameters(p);
 782     a.ResizeTo(4,4);
 783     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
 784     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
 785     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
 786     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 787     return true;
 788   }
 789 }
 790
 791 void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
 792                                         const AliExternalTrackParam &tp2,
 793                                         Int_t s1,Int_t s2) {
 794   //
 795   // Fill residual histograms
 796   // Innner-Outer
 797   // Left right - x-y
 798   // A-C side
 799   if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0)  {
 800     GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
 801     GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
 802     GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
 803     GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
 804   }
 805 }
 806
 807
 808
 809 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
 810 {
 811   //
 812   // return specified residual histogram - it is only QA
 813   // if force specified the histogram and given histogram is not existing
 814   //  new histogram is created
 815   //
 816   if (GetIndex(s1,s2)>=72*72) return 0;
 817   TObjArray *histoArray=0;
 818   switch (type) {
 819   case kY:
 820     histoArray = &fDyHistArray; break;
 821   case kZ:
 822     histoArray = &fDzHistArray; break;
 823   case kPhi:
 824     histoArray = &fDphiHistArray; break;
 825   case kTheta:
 826     histoArray = &fDthetaHistArray; break;
 827   }
 828   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
 829   if (histo) return histo;
 830   if (force==kFALSE) return 0;
 831   //
 832   stringstream name;
 833   stringstream title;
 834   switch (type) {
 835   case kY:
 836     name<<"hist_y_"<<s1<<"_"<<s2;
 837     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
 838     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
 839     break;
 840   case kZ:
 841     name<<"hist_z_"<<s1<<"_"<<s2;
 842     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
 843     histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
 844     break;
 845   case kPhi:
 846     name<<"hist_phi_"<<s1<<"_"<<s2;
 847     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
 848     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
 849     break;
 850   case kTheta:
 851     name<<"hist_theta_"<<s1<<"_"<<s2;
 852     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
 853     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
 854     break;
 855   }
 856   histo->SetDirectory(0);
 857   histoArray->AddAt(histo,GetIndex(s1,s2));
 858   return histo;
 859 }
 860
 861 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
 862                                            Int_t i0, Int_t i1, FitType type)
 863 {
 864   //
 865   //
 866   //
 867   TMatrixD mat;
 868   TObjArray *fitArray=0;
 869   Double_t xsec[1000];
 870   Double_t ysec[1000];
 871   Int_t npoints=0;
 872   for (Int_t isec = sec0; isec<=sec1; isec++){
 873     Int_t isec2 = (isec+dsec)%72;
 874     switch (type) {
 875     case k6:
 876       GetTransformation6(isec,isec2,mat);break;
 877     case k9:
 878       GetTransformation9(isec,isec2,mat);break;
 879     case k12:
 880       GetTransformation12(isec,isec2,mat);break;
 881     }
 882     xsec[npoints]=isec;
 883     ysec[npoints]=mat(i0,i1);
 884     ++npoints;
 885   }
 886   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
 887   Char_t name[1000];
 888   sprintf(name,"Mat[%d,%d]  Type=%d",i0,i1,type);
 889   gr->SetName(name);
 890   return gr;
 891 }
 892
 893 void  AliTPCcalibAlign::MakeTree(const char *fname){
 894   //
 895   // make tree with alignment cosntant  -
 896   // For  QA visualization
 897   //
 898   const Int_t kMinPoints=50;
 899   TTreeSRedirector cstream(fname);
 900   for (Int_t s1=0;s1<72;++s1)
 901     for (Int_t s2=0;s2<72;++s2){
 902       if (fPoints[GetIndex(s1,s2)]<kMinPoints) continue;
 903       TMatrixD m6;
 904       TMatrixD m9;
 905       TMatrixD m12;
 906       GetTransformation6(s1,s2,m6);
 907       GetTransformation9(s1,s2,m9);
 908       GetTransformation12(s1,s2,m12);
 909       Double_t dy=0, dz=0, dphi=0,dtheta=0;
 910       Double_t sy=0, sz=0, sphi=0,stheta=0;
 911       Double_t ny=0, nz=0, nphi=0,ntheta=0;
 912       TH1 * his=0;
 913       his = GetHisto(kY,s1,s2);
 914       if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
 915       his = GetHisto(kZ,s1,s2);
 916       if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
 917       his = GetHisto(kPhi,s1,s2);
 918       if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
 919       his = GetHisto(kTheta,s1,s2);
 920       if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
 921       //
 922       cstream<<"Align"<<
 923         "s1="<<s1<<     // reference sector
 924         "s2="<<s2<<     // sector to align
 925         "m6.="<<&m6<<   // tranformation matrix
 926         "m9.="<<&m9<<   //
 927         "m12.="<<&m12<<
 928         //               histograms mean RMS and entries
 929         "dy="<<dy<<
 930         "sy="<<sy<<
 931         "ny="<<ny<<
 932         "dz="<<dz<<
 933         "sz="<<sz<<
 934         "nz="<<nz<<
 935         "dphi="<<dphi<<
 936         "sphi="<<sphi<<
 937         "nphi="<<nphi<<
 938         "dtheta="<<dtheta<<
 939         "stheta="<<stheta<<
 940         "ntheta="<<ntheta<<
 941         "\n";
 942     }
 943
 944 }