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