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