TPC/AliTPCcalibAlign.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //     Class to make a internal alignemnt of TPC chambers                    //
  20 //
  21 //     Requierements - Warnings:
  22 //     1. Before using this componenent the magnetic filed has to be set properly //
  23 //     2. The systematic effects  - unlinearities has to be understood
  24 //
  25 //     If systematic and unlinearities are not under control
  26 //     the alignment is just effective alignment. Not second order corrction
  27 //     are calculated.
  28 //
  29 //     The histograming of the edge effects and unlineratities integral part
  30 //     of the component (currently only in debug stream)
  31 //
  32 //     3 general type of linear transformation investigated (see bellow)
  33 //
  34 //     By default only 6 parameter alignment to be used - other just for QA purposes
  35
  36 //     Different linear tranformation investigated
  37 //     12 parameters - arbitrary linear transformation
  38 //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
  39 //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
  40 //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
  41 //
  42 //      9 parameters - scaling fixed to 1
  43 //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
  44 //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
  45 //                     a20  a21  a22 a23     p[4]   p[5]  1      p[8]
  46 //
  47 //      6 parameters - x-y rotation x-z, y-z tiliting
  48 //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
  49 //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
  50 //                     a20  a21  a22 a23     p[1]   p[2]  1     p[5]
  51 //
  52 //
  53 //      Debug stream supported
  54 //      0. Align    - The main output of the Alignment component
  55 //                  - Used for visualization of the misalignment between sectors
  56 //                  - Results of the missalignment fit and the mean and sigmas of histograms
  57 //                   stored there
  58 //      1. Tracklet - StreamLevel >1
  59 //                  - Dump all information about tracklet match from sector1 to sector 2
  60 //                  - Default histogram residulas created in parallel
  61 //                  - Check this streamer in case of suspicious content of these histograms
  62 //      2. Track    - StreamLevel>5
  63 //                  - For debugging of the edge effects
  64 //                  - All information  - extrapolation inside of one sectors
  65 //                  - Created in order to distinguish between unlinearities inside of o
  66 //                    sector and  missalignment
  67
  68 //
  69 //
  70 /*
  71   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  72   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  73   AliXRDPROOFtoolkit tool;
  74   TChain * chain = tool.MakeChain("align.txt","Track",0,10200);
  75   chain->Lookup();
  76   TCut cutA("abs(tp1.fP[1]-tp2.fP[1])<0.3&&abs(tp1.fP[0]-tp2.fP[0])<0.15&&abs(tp1.fP[3]-tp2.fP[3])<0.01&&abs(tp1.fP[2]-tp2.fP[2])<0.01");
  77   TCut cutS("s1%36==s2%36");
  78
  79   .x ~/UliStyle.C
  80   .x $ALICE_ROOT/macros/loadlibsREC.C
  81
  82   gSystem->Load("$ROOTSYS/lib/libXrdClient.so");
  83   gSystem->Load("libProof");
  84   gSystem->Load("libANALYSIS");
  85   gSystem->Load("libSTAT");
  86   gSystem->Load("libTPCcalib");
  87   //
  88   // compare reference
  89   TFile fcalib("CalibObjects.root");
  90   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  91
  92   AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
  93   //
  94   //
  95   align->EvalFitters();
  96   align->MakeTree("alignTree.root");
  97   TFile falignTree("alignTree.root");
  98   TTree * treeAlign = (TTree*)falignTree.Get("Align");
  99
 100
 101 */
 102
 103 ////
 104 ////
 105
 106 #include "TLinearFitter.h"
 107 #include "AliTPCcalibAlign.h"
 108 #include "AliExternalTrackParam.h"
 109 #include "AliTPCTracklet.h"
 110 #include "TH1D.h"
 111 #include "TH2F.h"
 112 #include "TVectorD.h"
 113 #include "TTreeStream.h"
 114 #include "TFile.h"
 115 #include "TF1.h"
 116 #include "TGraphErrors.h"
 117 #include "AliTPCclusterMI.h"
 118 #include "AliTPCseed.h"
 119 #include "AliTracker.h"
 120 #include "TClonesArray.h"
 121
 122
 123 #include "TTreeStream.h"
 124 #include <iostream>
 125 #include <sstream>
 126 using namespace std;
 127
 128 ClassImp(AliTPCcalibAlign)
 129
 130 AliTPCcalibAlign::AliTPCcalibAlign()
 131   :  AliTPCcalibBase(),
 132      fDphiHistArray(72*72),
 133      fDthetaHistArray(72*72),
 134      fDyHistArray(72*72),
 135      fDzHistArray(72*72),
 136      //
 137      fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
 138      fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
 139      fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
 140      fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
 141      fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
 142      fDzZHistArray(72*72),        // array of residual histograms  z     -kZz
 143      fFitterArray12(72*72),
 144      fFitterArray9(72*72),
 145      fFitterArray6(72*72)
 146 {
 147   //
 148   // Constructor
 149   //
 150   for (Int_t i=0;i<72*72;++i) {
 151     fPoints[i]=0;
 152   }
 153 }
 154
 155 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
 156   :AliTPCcalibBase(),
 157    fDphiHistArray(72*72),
 158    fDthetaHistArray(72*72),
 159    fDyHistArray(72*72),
 160    fDzHistArray(72*72),
 161    fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
 162    fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
 163    fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
 164    fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
 165    fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
 166    fDzZHistArray(72*72),        // array of residual histograms  z     -kZz
 167
 168    fFitterArray12(72*72),
 169    fFitterArray9(72*72),
 170    fFitterArray6(72*72)
 171 {
 172   //
 173   // Constructor
 174   //
 175   SetName(name);
 176   SetTitle(title);
 177   for (Int_t i=0;i<72*72;++i) {
 178     fPoints[i]=0;
 179   }
 180 }
 181
 182
 183 AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
 184   :AliTPCcalibBase(align),
 185    fDphiHistArray(align.fDphiHistArray),
 186    fDthetaHistArray(align.fDthetaHistArray),
 187    fDyHistArray(align.fDyHistArray),
 188    fDzHistArray(align.fDzHistArray),
 189    fDyPhiHistArray(align.fDyPhiHistArray),      // array of residual histograms  y     -kYPhi
 190    fDzThetaHistArray(align.fDzThetaHistArray),    // array of residual histograms  z-z   -kZTheta
 191    fDphiZHistArray(align.fDphiZHistArray),      // array of residual histograms  phi   -kPhiz
 192    fDthetaZHistArray(align.fDthetaZHistArray),    // array of residual histograms  theta -kThetaz
 193    fDyZHistArray(align.fDyZHistArray),        // array of residual histograms  y     -kYz
 194    fDzZHistArray(align.fDzZHistArray),        // array of residual histograms  z     -kZz
 195    //
 196    fFitterArray12(align.fFitterArray12),
 197    fFitterArray9(align.fFitterArray9),
 198    fFitterArray6(align.fFitterArray6)
 199
 200 {
 201   //
 202   // copy constructor - copy also the content
 203   //
 204   TH1 * his = 0;
 205   TObjArray * arr0=0;
 206   const TObjArray *arr1=0;
 207   for (Int_t index =0; index<72*72; index++){
 208     for (Int_t iarray=0;iarray<10; iarray++){
 209       if (iarray==kY){
 210         arr0 = &fDyHistArray;
 211         arr1 = &align.fDyHistArray;
 212       }
 213       if (iarray==kZ){
 214         arr0 = &fDzHistArray;
 215         arr1 = &align.fDzHistArray;
 216       }
 217       if (iarray==kPhi){
 218         arr0 = &fDphiHistArray;
 219         arr1 = &align.fDphiHistArray;
 220       }
 221       if (iarray==kTheta){
 222         arr0 = &fDthetaHistArray;
 223         arr1 = &align.fDthetaHistArray;
 224       }
 225       if (iarray==kYz){
 226         arr0 = &fDyZHistArray;
 227         arr1 = &align.fDyZHistArray;
 228       }
 229       if (iarray==kZz){
 230         arr0 = &fDzZHistArray;
 231         arr1 = &align.fDzZHistArray;
 232       }
 233       if (iarray==kPhiZ){
 234         arr0 = &fDphiZHistArray;
 235         arr1 = &align.fDphiZHistArray;
 236       }
 237       if (iarray==kThetaZ){
 238         arr0 = &fDthetaZHistArray;
 239         arr1 = &align.fDthetaZHistArray;
 240       }
 241
 242       if (iarray==kYPhi){
 243         arr0 = &fDyPhiHistArray;
 244         arr1 = &align.fDyPhiHistArray;
 245       }
 246       if (iarray==kZTheta){
 247         arr0 = &fDzThetaHistArray;
 248         arr1 = &align.fDzThetaHistArray;
 249       }
 250
 251       if (arr1->At(index)) {
 252         his = (TH1*)arr1->At(index)->Clone();
 253         his->SetDirectory(0);
 254         arr0->AddAt(his,index);
 255       }
 256     }
 257   }
 258 }
 259
 260
 261 AliTPCcalibAlign::~AliTPCcalibAlign() {
 262   //
 263   // destructor
 264   //
 265 }
 266
 267 void AliTPCcalibAlign::Process(AliTPCseed *seed) {
 268   //
 269   //
 270   //
 271
 272   TObjArray tracklets=
 273     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
 274                                     kFALSE,20,2);
 275   //  TObjArray trackletsL=
 276 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
 277 //                                  kFALSE,20,2);
 278 //   TObjArray trackletsQ=
 279 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
 280 //                                  kFALSE,20,2);
 281 //   TObjArray trackletsR=
 282 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
 283 //                                  kFALSE,20,2);
 284   tracklets.SetOwner();
 285  //  trackletsL.SetOwner();
 286 //   trackletsQ.SetOwner();
 287 //   trackletsR.SetOwner();
 288   if (tracklets.GetEntries()==2) {
 289     AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
 290     AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
 291     if (t1->GetSector()>t2->GetSector()) {
 292       AliTPCTracklet* tmp=t1;
 293       t1=t2;
 294       t2=tmp;
 295     }
 296     AliExternalTrackParam *common1=0,*common2=0;
 297     if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
 298       ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
 299     delete common1;
 300     delete common2;
 301   }
 302
 303 }
 304
 305 void AliTPCcalibAlign::Analyze(){
 306   //
 307   // Analyze function
 308   //
 309   EvalFitters();
 310 }
 311
 312
 313 void AliTPCcalibAlign::Terminate(){
 314   //
 315   // Terminate function
 316   // call base terminate + Eval of fitters
 317   //
 318   Info("AliTPCcalibAlign","Terminate");
 319   EvalFitters();
 320   AliTPCcalibBase::Terminate();
 321 }
 322
 323
 324
 325
 326 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
 327                                         const AliExternalTrackParam &tp2,
 328                                         const AliTPCseed * seed,
 329                                         Int_t s1,Int_t s2) {
 330
 331   //
 332   //
 333   //
 334   //
 335   //
 336   // Process function to fill fitters
 337   //
 338   Double_t t1[5],t2[5];
 339   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 340   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 341   x1   =tp1.GetX();
 342   y1   =tp1.GetY();
 343   z1   =tp1.GetZ();
 344   Double_t snp1=tp1.GetSnp();
 345   dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
 346   Double_t tgl1=tp1.GetTgl();
 347   // dz/dx = 1/(cos(theta)*cos(phi))
 348   dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1);
 349   x2   =tp2.GetX();
 350   y2   =tp2.GetY();
 351   z2   =tp2.GetZ();
 352   Double_t snp2=tp2.GetSnp();
 353   dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
 354   Double_t tgl2=tp2.GetTgl();
 355   dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2);
 356   //
 357   //
 358   //
 359   if (fStreamLevel>1){
 360     TTreeSRedirector *cstream = GetDebugStreamer();
 361     if (cstream){
 362       static TVectorD vec1(5);
 363       static TVectorD vec2(5);
 364       vec1.SetElements(t1);
 365       vec2.SetElements(t2);
 366       AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
 367       AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
 368       (*cstream)<<"Tracklet"<<
 369         "run="<<fRun<<              //  run number
 370         "event="<<fEvent<<          //  event number
 371         "time="<<fTime<<            //  time stamp of event
 372         "trigger="<<fTrigger<<      //  trigger
 373         "mag="<<fMagF<<             //  magnetic field
 374         "tp1.="<<p1<<
 375         "tp2.="<<p2<<
 376         "v1.="<<&vec1<<
 377         "v2.="<<&vec2<<
 378         "s1="<<s1<<
 379         "s2="<<s2<<
 380         "\n";
 381     }
 382   }
 383   //
 384   // Aplly cut selection
 385   /*
 386     TChain * chainalign = tool.MakeChain("align.txt","Tracklet",0,1000000)
 387     chainalign->Lookup();
 388
 389     // Cuts to be justified with the debug streamer
 390     //
 391     TCut c1pt("abs((tp1.fP[4]+tp2.fP[4])*0.5)<3"); // pt cut  - OK
 392     TCut cdy("abs(tp1.fP[0]-tp2.fP[0])<2");
 393     TCut cdz("abs(tp1.fP[1]-tp2.fP[1])<2");
 394     TCut cdphi("abs(tp1.fP[2]-tp2.fP[2])<0.02");
 395     TCut cdt("abs(tp1.fP[3]-tp2.fP[3])<0.02");
 396     TCut cd1pt("abs(tp1.fP[4]-tp2.fP[4])<0.3");    // delta 1/pt cut  -OK
 397     //
 398     //
 399     TCut acut =  c1pt+cdy+cdz+cdphi+cdt+cd1pt;
 400     chainalign->Draw(">>listEL",acut,"entryList");
 401     TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
 402     chainalign->SetEntryList(elist);
 403
 404   */
 405   //   1. pt cut
 406   //   2. dy
 407   //   3. dz
 408   //   4. dphi
 409   //   5. dtheta
 410   //   6. d1pt
 411   if (GetDebugLevel()>50) printf("Process track\n");
 412   if (TMath::Abs(tp1.GetParameter()[0]-tp2.GetParameter()[0])>2)    return;
 413   if (TMath::Abs(tp1.GetParameter()[1]-tp2.GetParameter()[1])>2)    return;
 414   if (TMath::Abs(tp1.GetParameter()[2]-tp2.GetParameter()[2])>0.02) return;
 415   if (TMath::Abs(tp1.GetParameter()[3]-tp2.GetParameter()[3])>0.02) return;
 416   if (TMath::Abs(tp1.GetParameter()[4]-tp2.GetParameter()[4])>0.3)  return;
 417   if (TMath::Abs((tp1.GetParameter()[4]+tp2.GetParameter()[4])*0.5)>3)  return;
 418   if (TMath::Abs((tp1.GetParameter()[0]-tp2.GetParameter()[0]))<0.000000001)  return;
 419    if (GetDebugLevel()>50) printf("Filling track\n");
 420   //
 421   // fill resolution histograms - previous cut included
 422   FillHisto(tp1,tp2,s1,s2);
 423   //
 424   Process12(t1,t2,GetOrMakeFitter12(s1,s2));
 425   Process9(t1,t2,GetOrMakeFitter9(s1,s2));
 426   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
 427   ProcessDiff(tp1,tp2, seed,s1,s2);
 428   ++fPoints[GetIndex(s1,s2)];
 429 }
 430
 431 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
 432                                     const AliExternalTrackParam &t2,
 433                                     const AliTPCseed *seed,
 434                                     Int_t s1,Int_t s2)
 435 {
 436   //
 437   // Process local residuals function
 438   //
 439   TVectorD vecX(160);
 440   TVectorD vecY(160);
 441   TVectorD vecZ(160);
 442   TVectorD vecClY(160);
 443   TVectorD vecClZ(160);
 444   TClonesArray arrCl("AliTPCclusterMI",160);
 445   arrCl.ExpandCreateFast(160);
 446   Int_t count1=0, count2=0;
 447   for (Int_t i=0;i<160;++i) {
 448     AliTPCclusterMI *c=seed->GetClusterPointer(i);
 449     vecX[i]=0;
 450     vecY[i]=0;
 451     vecZ[i]=0;
 452     if (!c) continue;
 453     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
 454     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
 455     vecClY[i] = c->GetY();
 456     vecClZ[i] = c->GetZ();
 457     cl=*c;
 458     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
 459     if (c->GetDetector()==s1) ++count1;
 460     if (c->GetDetector()==s2) ++count2;
 461     Double_t gxyz[3],xyz[3];
 462     t1.GetXYZ(gxyz);
 463     Float_t bz = AliTracker::GetBz(gxyz);
 464     par->GetYAt(c->GetX(), bz, xyz[1]);
 465     par->GetZAt(c->GetX(), bz, xyz[2]);
 466     vecX[i] = c->GetX();
 467     vecY[i]= xyz[1];
 468     vecZ[i]= xyz[2];
 469   }
 470   //
 471   //
 472   if (fStreamLevel>5){
 473     //
 474     // huge output - cluster residuals to be investigated
 475     //
 476     TTreeSRedirector *cstream = GetDebugStreamer();
 477     //AliTPCseed * t = (AliTPCseed*) seed;
 478     //AliExternalTrackParam *p0 = &((AliExternalTrackParam&)seed);
 479     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
 480     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
 481     /*
 482
 483       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");
 484
 485     */
 486
 487     if (cstream){
 488       (*cstream)<<"Track"<<
 489         "run="<<fRun<<              //  run number
 490         "event="<<fEvent<<          //  event number
 491         "time="<<fTime<<            //  time stamp of event
 492         "trigger="<<fTrigger<<      //  trigger
 493         "mag="<<fMagF<<             //  magnetic field
 494         "Cl.="<<&arrCl<<
 495         //"tp0.="<<p0<<
 496         "tp1.="<<p1<<
 497         "tp2.="<<p2<<
 498         "vtX.="<<&vecX<<
 499         "vtY.="<<&vecY<<
 500         "vtZ.="<<&vecZ<<
 501         "vcY.="<<&vecClY<<
 502         "vcZ.="<<&vecClZ<<
 503         "s1="<<s1<<
 504         "s2="<<s2<<
 505         "c1="<<count1<<
 506         "c2="<<count2<<
 507         "\n";
 508     }
 509   }
 510 }
 511
 512
 513
 514
 515 void AliTPCcalibAlign::Process12(const Double_t *t1,
 516                                  const Double_t *t2,
 517                                  TLinearFitter *fitter) {
 518   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 519   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 520   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 521   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 522   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 523   //
 524   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
 525   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
 526   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11]
 527
 528
 529
 530   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 531   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 532
 533   // TODO:
 534   Double_t sy    = 0.1;
 535   Double_t sz    = 0.1;
 536   Double_t sdydx = 0.001;
 537   Double_t sdzdx = 0.001;
 538
 539   Double_t p[12];
 540   Double_t value;
 541
 542   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 543   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 544   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 545   for (Int_t i=0; i<12;i++) p[i]=0.;
 546   p[3+0] = x1;          // a10
 547   p[3+1] = y1;          // a11
 548   p[3+2] = z1;          // a12
 549   p[9+1] = 1.;          // a13
 550   p[0+1] = y1*dydx2;    // a01
 551   p[0+2] = z1*dydx2;    // a02
 552   p[9+0] = dydx2;       // a03
 553   value  = y2;
 554   fitter->AddPoint(p,value,sy);
 555
 556   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 557   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 558   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 559   for (Int_t i=0; i<12;i++) p[i]=0.;
 560   p[6+0] = x1;           // a20
 561   p[6+1] = y1;           // a21
 562   p[6+2] = z1;           // a22
 563   p[9+2] = 1.;           // a23
 564   p[0+1] = y1*dzdx2;     // a01
 565   p[0+2] = z1*dzdx2;     // a02
 566   p[9+0] = dzdx2;        // a03
 567   value  = z2;
 568   fitter->AddPoint(p,value,sz);
 569
 570   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 571   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 572   for (Int_t i=0; i<12;i++) p[i]=0.;
 573   p[3+0] = 1.;           // a10
 574   p[3+1] = dydx1;        // a11
 575   p[3+2] = dzdx1;        // a12
 576   p[0+0] = -dydx2;       // a00
 577   p[0+1] = -dydx1*dydx2; // a01
 578   p[0+2] = -dzdx1*dydx2; // a02
 579   value  = 0.;
 580   fitter->AddPoint(p,value,sdydx);
 581
 582   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 583   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 584   for (Int_t i=0; i<12;i++) p[i]=0.;
 585   p[6+0] = 1;            // a20
 586   p[6+1] = dydx1;        // a21
 587   p[6+2] = dzdx1;        // a22
 588   p[0+0] = -dzdx2;       // a00
 589   p[0+1] = -dydx1*dzdx2; // a01
 590   p[0+2] = -dzdx1*dzdx2; // a02
 591   value  = 0.;
 592   fitter->AddPoint(p,value,sdzdx);
 593 }
 594
 595 void AliTPCcalibAlign::Process9(Double_t *t1,
 596                                 Double_t *t2,
 597                                 TLinearFitter *fitter) {
 598   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
 599   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
 600   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
 601   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 602   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 603   //
 604   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
 605   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
 606   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8]
 607
 608
 609   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 610   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 611
 612   // TODO:
 613   Double_t sy    = 0.1;
 614   Double_t sz    = 0.1;
 615   Double_t sdydx = 0.001;
 616   Double_t sdzdx = 0.001;
 617   //
 618   Double_t p[12];
 619   Double_t value;
 620
 621   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 622   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 623   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 624   for (Int_t i=0; i<12;i++) p[i]=0.;
 625   p[2]   += x1;           // a10
 626   //p[]  +=1;             // a11
 627   p[3]   += z1;           // a12
 628   p[7]   += 1;            // a13
 629   p[0]   += y1*dydx2;     // a01
 630   p[1]   += z1*dydx2;     // a02
 631   p[6]   += dydx2;        // a03
 632   value   = y2-y1;        //-a11
 633   fitter->AddPoint(p,value,sy);
 634   //
 635   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 636   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 637   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 638   for (Int_t i=0; i<12;i++) p[i]=0.;
 639   p[4]   += x1;           // a20
 640   p[5]   += y1;           // a21
 641   //p[]  += z1;           // a22
 642   p[8]   += 1.;           // a23
 643   p[0]   += y1*dzdx2;     // a01
 644   p[1]   += z1*dzdx2;     // a02
 645   p[6]   += dzdx2;        // a03
 646   value  = z2-z1;         //-a22
 647   fitter->AddPoint(p,value,sz);
 648
 649   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 650   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 651   for (Int_t i=0; i<12;i++) p[i]=0.;
 652   p[2]   += 1.;           // a10
 653   //p[]  += dydx1;      // a11
 654   p[3]   += dzdx1;        // a12
 655   //p[]  += -dydx2;       // a00
 656   p[0]   += -dydx1*dydx2; // a01
 657   p[1]   += -dzdx1*dydx2; // a02
 658   value  = -dydx1+dydx2;  // -a11 + a00
 659   fitter->AddPoint(p,value,sdydx);
 660
 661   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 662   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 663   for (Int_t i=0; i<12;i++) p[i]=0.;
 664   p[4]   += 1;            // a20
 665   p[5]   += dydx1;        // a21
 666   //p[]  += dzdx1;        // a22
 667   //p[]  += -dzdx2;       // a00
 668   p[0]   += -dydx1*dzdx2; // a01
 669   p[1]   += -dzdx1*dzdx2; // a02
 670   value  = -dzdx1+dzdx2;  // -a22 + a00
 671   fitter->AddPoint(p,value,sdzdx);
 672 }
 673
 674 void AliTPCcalibAlign::Process6(Double_t *t1,
 675                                 Double_t *t2,
 676                                 TLinearFitter *fitter) {
 677   // x2    =  1  *x1 +-a01*y1 + 0      +a03
 678   // y2    =  a01*x1 + 1  *y1 + 0      +a13
 679   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
 680   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 681   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
 682   //
 683   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
 684   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
 685   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5]
 686
 687   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
 688   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 689
 690   // TODO:
 691   Double_t sy    = 0.1;
 692   Double_t sz    = 0.1;
 693   Double_t sdydx = 0.001;
 694   Double_t sdzdx = 0.001;
 695
 696   Double_t p[12];
 697   Double_t value;
 698   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 699   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
 700   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
 701   for (Int_t i=0; i<12;i++) p[i]=0.;
 702   p[0]   += x1;           // a10
 703   //p[]  +=1;             // a11
 704   //p[]  += z1;           // a12
 705   p[4]   += 1;            // a13
 706   p[0]   += -y1*dydx2;    // a01
 707   //p[]  += z1*dydx2;     // a02
 708   p[3]   += dydx2;        // a03
 709   value   = y2-y1;        //-a11
 710   fitter->AddPoint(p,value,sy);
 711   //
 712   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
 713   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
 714   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
 715   for (Int_t i=0; i<12;i++) p[i]=0.;
 716   p[1]   += x1;           // a20
 717   p[2]   += y1;           // a21
 718   //p[]  += z1;           // a22
 719   p[5]   += 1.;           // a23
 720   p[0]   += -y1*dzdx2;    // a01
 721   //p[]   += z1*dzdx2;     // a02
 722   p[3]   += dzdx2;        // a03
 723   value  = z2-z1;         //-a22
 724   fitter->AddPoint(p,value,sz);
 725
 726   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 727   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
 728   for (Int_t i=0; i<12;i++) p[i]=0.;
 729   p[0]   += 1.;           // a10
 730   //p[]  += dydx1;      // a11
 731   //p[]   += dzdx1;        // a12
 732   //p[]  += -dydx2;       // a00
 733   p[0]   +=  dydx1*dydx2; // a01
 734   //p[]   += -dzdx1*dydx2; // a02
 735   value  = -dydx1+dydx2;  // -a11 + a00
 736   fitter->AddPoint(p,value,sdydx);
 737
 738   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
 739   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
 740   for (Int_t i=0; i<12;i++) p[i]=0.;
 741   p[1]   += 1;            // a20
 742   p[2]   += dydx1;        // a21
 743   //p[]  += dzdx1;        // a22
 744   //p[]  += -dzdx2;       // a00
 745   p[0]   +=  dydx1*dzdx2; // a01
 746   //p[]  += -dzdx1*dzdx2; // a02
 747   value  = -dzdx1+dzdx2;  // -a22 + a00
 748   fitter->AddPoint(p,value,sdzdx);
 749 }
 750
 751
 752
 753
 754 void AliTPCcalibAlign::EvalFitters() {
 755   //
 756   // Analyze function
 757   //
 758   // Perform the fitting using linear fitters
 759   //
 760   Int_t kMinPoints =50;
 761   TLinearFitter *f;
 762   TFile fff("alignDebug.root","recreate");
 763   for (Int_t s1=0;s1<72;++s1)
 764     for (Int_t s2=0;s2<72;++s2){
 765       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 766         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 767         if (f->Eval()!=0) {
 768           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 769           f->Write(Form("f12_%d_%d",s1,s2));
 770         }else{
 771           f->Write(Form("f12_%d_%d",s1,s2));
 772         }
 773       }
 774       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 775         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 776         if (f->Eval()!=0) {
 777           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 778         }else{
 779           f->Write(Form("f9_%d_%d",s1,s2));
 780         }
 781       }
 782       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
 783         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
 784         if (f->Eval()!=0) {
 785           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
 786         }else{
 787           f->Write(Form("f6_%d_%d",s1,s2));
 788         }
 789       }
 790     }
 791   this->Write("align");
 792   /*
 793
 794   fitter->Eval();
 795   fitter->Eval();
 796   chi212 = align->GetChisquare()/(4.*entries);
 797
 798   TMatrixD mat(13,13);
 799   TVectorD par(13);
 800   align->GetParameters(par);
 801   align->GetCovarianceMatrix(mat);
 802
 803   //
 804   //
 805   for (Int_t i=0; i<12;i++){
 806     palign12(i)= par(i+1);
 807     for (Int_t j=0; j<12;j++){
 808       pcovar12(i,j)   = mat(i+1,j+1);
 809       pcovar12(i,j) *= chi212;
 810     }
 811   }
 812   //
 813   for (Int_t i=0; i<12;i++){
 814     psigma12(i)  = TMath::Sqrt(pcovar12(i,i));
 815     palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i));
 816     for (Int_t j=0; j<12;j++){
 817       pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j));
 818     }
 819   }
 820   */
 821 }
 822
 823 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
 824   //
 825   // get or make fitter - general linear transformation
 826   //
 827   static Int_t counter12=0;
 828   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]");
 829   TLinearFitter * fitter = GetFitter12(s1,s2);
 830   if (fitter) return fitter;
 831   //  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]");
 832   fitter =new TLinearFitter(&f12,"");
 833   fitter->StoreData(kTRUE);
 834   fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
 835   counter12++;
 836   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
 837   return fitter;
 838 }
 839
 840 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
 841   //
 842   //get or make fitter - general linear transformation - no scaling
 843   //
 844   static Int_t counter9=0;
 845   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 846   TLinearFitter * fitter = GetFitter9(s1,s2);
 847   if (fitter) return fitter;
 848   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
 849   fitter =new TLinearFitter(&f9,"");
 850   fitter->StoreData(kTRUE);
 851   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
 852   counter9++;
 853   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
 854   return fitter;
 855 }
 856
 857 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
 858   //
 859   // get or make fitter  - 6 paramater linear tranformation
 860   //                     - no scaling
 861   //                     - rotation x-y
 862   //                     - tilting x-z, y-z
 863   static Int_t counter6=0;
 864   static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
 865   TLinearFitter * fitter = GetFitter6(s1,s2);
 866   if (fitter) return fitter;
 867   //  fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
 868   fitter=new TLinearFitter(&f6,"");
 869   fitter->StoreData(kTRUE);
 870   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
 871   counter6++;
 872   if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
 873   return fitter;
 874 }
 875
 876
 877
 878
 879
 880 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
 881   //
 882   // GetTransformation matrix - 12 paramaters - generael linear transformation
 883   //
 884   if (!GetFitter12(s1,s2))
 885     return false;
 886   else {
 887     TVectorD p(12);
 888     GetFitter12(s1,s2)->GetParameters(p);
 889     a.ResizeTo(4,4);
 890     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
 891     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
 892     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
 893     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 894     return true;
 895   }
 896 }
 897
 898 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
 899   //
 900   // GetTransformation matrix - 9 paramaters - general linear transformation
 901   //                            No scaling
 902   //
 903   if (!GetFitter9(s1,s2))
 904     return false;
 905   else {
 906     TVectorD p(9);
 907     GetFitter9(s1,s2)->GetParameters(p);
 908     a.ResizeTo(4,4);
 909     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
 910     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
 911     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
 912     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 913     return true;
 914   }
 915 }
 916
 917 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
 918   //
 919   // GetTransformation matrix - 6  paramaters
 920   //                            3  translation
 921   //                            1  rotation -x-y
 922   //                            2  tilting x-z y-z
 923   if (!GetFitter6(s1,s2))
 924     return false;
 925   else {
 926     TVectorD p(6);
 927     GetFitter6(s1,s2)->GetParameters(p);
 928     a.ResizeTo(4,4);
 929     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
 930     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
 931     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
 932     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
 933     return true;
 934   }
 935 }
 936
 937 void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
 938                                         const AliExternalTrackParam &tp2,
 939                                         Int_t s1,Int_t s2) {
 940   //
 941   // Fill residual histograms
 942   // Innner-Outer
 943   // Left right - x-y
 944   // A-C side
 945   if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0)  {
 946     GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
 947     GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
 948     GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
 949     GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
 950     //
 951     GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
 952     GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
 953     GetHisto(kYz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetY()-tp2.GetY());
 954     GetHisto(kZz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetZ()-tp2.GetZ());
 955     //
 956     GetHisto(kYPhi,s1,s2,kTRUE)->Fill(tp1.GetSnp(),tp1.GetY()-tp2.GetY());
 957     GetHisto(kZTheta,s1,s2,kTRUE)->Fill(tp1.GetTgl(),tp1.GetZ()-tp2.GetZ());
 958
 959
 960   }
 961 }
 962
 963
 964
 965 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
 966 {
 967   //
 968   // return specified residual histogram - it is only QA
 969   // if force specified the histogram and given histogram is not existing
 970   //  new histogram is created
 971   //
 972   if (GetIndex(s1,s2)>=72*72) return 0;
 973   TObjArray *histoArray=0;
 974   switch (type) {
 975   case kY:
 976     histoArray = &fDyHistArray; break;
 977   case kZ:
 978     histoArray = &fDzHistArray; break;
 979   case kPhi:
 980     histoArray = &fDphiHistArray; break;
 981   case kTheta:
 982     histoArray = &fDthetaHistArray; break;
 983   case kYPhi:
 984     histoArray = &fDyPhiHistArray; break;
 985   case kZTheta:
 986     histoArray = &fDzThetaHistArray; break;
 987   case kYz:
 988     histoArray = &fDyZHistArray; break;
 989   case kZz:
 990     histoArray = &fDzZHistArray; break;
 991   case kPhiZ:
 992     histoArray = &fDphiZHistArray; break;
 993   case kThetaZ:
 994     histoArray = &fDthetaZHistArray; break;
 995   }
 996   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
 997   if (histo) return histo;
 998   if (force==kFALSE) return 0;
 999   //
1000   stringstream name;
1001   stringstream title;
1002   switch (type) {
1003   case kY:
1004     name<<"hist_y_"<<s1<<"_"<<s2;
1005     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1006     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1007     break;
1008   case kZ:
1009     name<<"hist_z_"<<s1<<"_"<<s2;
1010     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1011     histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1012     break;
1013   case kPhi:
1014     name<<"hist_phi_"<<s1<<"_"<<s2;
1015     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1016     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1017     break;
1018   case kTheta:
1019     name<<"hist_theta_"<<s1<<"_"<<s2;
1020     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1021     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1022     break;
1023     //
1024     //
1025   case kYPhi:
1026     name<<"hist_yphi_"<<s1<<"_"<<s2;
1027     title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1028     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,128,-0.3,0.3); // +/- 3 mm
1029     break;
1030   case kZTheta:
1031     name<<"hist_ztheta_"<<s1<<"_"<<s2;
1032     title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1033     histo = new TH2F(name.str().c_str(),title.str().c_str(),128,20,-1,1,-0.3,0.3); // +/- 3 mm
1034     break;
1035     //
1036     //
1037     //
1038   case kYz:
1039     name<<"hist_yz_"<<s1<<"_"<<s2;
1040     title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1041     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1042     break;
1043   case kZz:
1044     name<<"hist_zz_"<<s1<<"_"<<s2;
1045     title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1046     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1047     break;
1048   case kPhiZ:
1049     name<<"hist_phiz_"<<s1<<"_"<<s2;
1050     title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1051     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1052     break;
1053   case kThetaZ:
1054     name<<"hist_thetaz_"<<s1<<"_"<<s2;
1055     title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1056     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1057     break;
1058
1059
1060   }
1061   histo->SetDirectory(0);
1062   histoArray->AddAt(histo,GetIndex(s1,s2));
1063   return histo;
1064 }
1065
1066 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
1067                                            Int_t i0, Int_t i1, FitType type)
1068 {
1069   //
1070   //
1071   //
1072   TMatrixD mat;
1073   //TObjArray *fitArray=0;
1074   Double_t xsec[1000];
1075   Double_t ysec[1000];
1076   Int_t npoints=0;
1077   for (Int_t isec = sec0; isec<=sec1; isec++){
1078     Int_t isec2 = (isec+dsec)%72;
1079     switch (type) {
1080     case k6:
1081       GetTransformation6(isec,isec2,mat);break;
1082     case k9:
1083       GetTransformation9(isec,isec2,mat);break;
1084     case k12:
1085       GetTransformation12(isec,isec2,mat);break;
1086     }
1087     xsec[npoints]=isec;
1088     ysec[npoints]=mat(i0,i1);
1089     ++npoints;
1090   }
1091   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1092   Char_t name[1000];
1093   sprintf(name,"Mat[%d,%d]  Type=%d",i0,i1,type);
1094   gr->SetName(name);
1095   return gr;
1096 }
1097
1098 void  AliTPCcalibAlign::MakeTree(const char *fname){
1099   //
1100   // make tree with alignment cosntant  -
1101   // For  QA visualization
1102   //
1103   /*
1104    TFile f("CalibObjects.root");
1105    TObjArray *array  = (TObjArray*)f.Get("TPCCalib");
1106    AliTPCcalibAlign   *alignTPC = (AliTPCcalibAlign   *)array->At(0);
1107    alignTPC->MakeTree("alignTree.root");
1108    TFile falign("alignTree.root");
1109    Align->Draw("dy")
1110    */
1111   const Int_t kMinPoints=50;
1112   TTreeSRedirector cstream(fname);
1113   for (Int_t s1=0;s1<72;++s1)
1114     for (Int_t s2=0;s2<72;++s2){
1115       if (fPoints[GetIndex(s1,s2)]<kMinPoints) continue;
1116       TMatrixD m6;
1117       TMatrixD m9;
1118       TMatrixD m12;
1119       GetTransformation6(s1,s2,m6);
1120       GetTransformation9(s1,s2,m9);
1121       GetTransformation12(s1,s2,m12);
1122       Double_t dy=0, dz=0, dphi=0,dtheta=0;
1123       Double_t sy=0, sz=0, sphi=0,stheta=0;
1124       Double_t ny=0, nz=0, nphi=0,ntheta=0;
1125       TH1 * his=0;
1126       his = GetHisto(kY,s1,s2);
1127       if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1128       his = GetHisto(kZ,s1,s2);
1129       if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1130       his = GetHisto(kPhi,s1,s2);
1131       if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
1132       his = GetHisto(kTheta,s1,s2);
1133       if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
1134       //
1135       cstream<<"Align"<<
1136         "s1="<<s1<<     // reference sector
1137         "s2="<<s2<<     // sector to align
1138         "m6.="<<&m6<<   // tranformation matrix
1139         "m9.="<<&m9<<   //
1140         "m12.="<<&m12<<
1141         //               histograms mean RMS and entries
1142         "dy="<<dy<<
1143         "sy="<<sy<<
1144         "ny="<<ny<<
1145         "dz="<<dz<<
1146         "sz="<<sz<<
1147         "nz="<<nz<<
1148         "dphi="<<dphi<<
1149         "sphi="<<sphi<<
1150         "nphi="<<nphi<<
1151         "dtheta="<<dtheta<<
1152         "stheta="<<stheta<<
1153         "ntheta="<<ntheta<<
1154         "\n";
1155     }
1156
1157 }
1158
1159
1160 //_____________________________________________________________________
1161 Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
1162   //
1163   // merge function
1164   //
1165   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
1166   if (!list)
1167     return 0;
1168   if (list->IsEmpty())
1169     return 1;
1170
1171   TIterator* iter = list->MakeIterator();
1172   TObject* obj = 0;
1173   iter->Reset();
1174   Int_t count=0;
1175   while((obj = iter->Next()) != 0)
1176     {
1177       AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
1178       if (entry == 0) continue;
1179       Add(entry);
1180       count++;
1181     }
1182   return count;
1183 }
1184
1185
1186 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
1187   //
1188   // Add entry - used for merging of compoents
1189   //
1190   for (Int_t i=0; i<72;i++){
1191     for (Int_t j=0; j<72;j++){
1192       if (align->fPoints[GetIndex(i,j)]==0) continue;
1193       fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
1194       //
1195       //
1196       //
1197       for (Int_t itype=0; itype<10; itype++){
1198         TH1 * his0=0, *his1=0;
1199         his0 = GetHisto((HistoType)itype,i,j);
1200         his1 = align->GetHisto((HistoType)itype,i,j);
1201         if (his1){
1202           if (his0) his0->Add(his1);
1203           else {
1204             his0 = GetHisto(kY,i,j,kTRUE);
1205             his0->Add(his1);
1206           }
1207         }
1208       }
1209     //   //
1210 //       // dy
1211 //       TH1* hdy0 = GetHisto(kY,i,j);
1212 //       TH1* hdy1 = align->GetHisto(kY,i,j);
1213 //       if (hdy1){
1214 //      if (hdy0) hdy0->Add(hdy1);
1215 //      else {
1216 //        hdy0 = GetHisto(kY,i,j,kTRUE);
1217 //        hdy0->Add(hdy1);
1218 //      }
1219 //       }
1220 //       //
1221 //       // dz
1222 //       TH1* hdz0 = GetHisto(kZ,i,j);
1223 //       TH1* hdz1 = align->GetHisto(kZ,i,j);
1224 //       if (hdz1){
1225 //      if (hdz0) hdz0->Add(hdz1);
1226 //      else {
1227 //        hdz0 = GetHisto(kZ,i,j,kTRUE);
1228 //        hdz0->Add(hdz1);
1229 //      }
1230 //       }
1231 //       //
1232 //       // dphi
1233 //       TH1* hdphi0 = GetHisto(kPhi,i,j);
1234 //       TH1* hdphi1 = align->GetHisto(kPhi,i,j);
1235 //       if (hdphi1){
1236 //      if (hdphi0) hdphi0->Add(hdphi1);
1237 //      else {
1238 //        hdphi0 = GetHisto(kPhi,i,j,kTRUE);
1239 //        hdphi0->Add(hdphi1);
1240 //      }
1241 //       }
1242 //       //
1243 //       // dtheta
1244 //       TH1* hdTheta0 = GetHisto(kTheta,i,j);
1245 //       TH1* hdTheta1 = align->GetHisto(kTheta,i,j);
1246 //       if (hdTheta1){
1247 //      if (hdTheta0) hdTheta0->Add(hdTheta1);
1248 //      else {
1249 //        hdTheta0 = GetHisto(kTheta,i,j,kTRUE);
1250 //        hdTheta0->Add(hdTheta1);
1251 //      }
1252 //       }
1253     }
1254   }
1255   TLinearFitter *f0=0;
1256   TLinearFitter *f1=0;
1257   for (Int_t i=0; i<72;i++){
1258     for (Int_t j=0; j<72;j++){
1259       if (align->fPoints[GetIndex(i,j)]==0) continue;
1260       //
1261       // fitter12
1262       f0 =  GetFitter12(i,j);
1263       f1 =  align->GetFitter12(i,j);
1264       if (f1){
1265         if (f0) f0->Add(f1);
1266         else {
1267           f0 = GetOrMakeFitter12(i,j);
1268           f0->Add(f1);
1269         }
1270       }
1271       //
1272       // fitter9
1273       f0 =  GetFitter9(i,j);
1274       f1 =  align->GetFitter9(i,j);
1275       if (f1){
1276         if (f0) f0->Add(f1);
1277         else {
1278           f0 = GetOrMakeFitter9(i,j);
1279           f0->Add(f1);
1280         }
1281       }
1282       f0 =  GetFitter6(i,j);
1283       f1 =  align->GetFitter6(i,j);
1284       if (f1){
1285         if (f0) f0->Add(f1);
1286         else {
1287           f0 = GetOrMakeFitter6(i,j);
1288           f0->Add(f1);
1289         }
1290       }
1291     }
1292   }
1293 }
1294
1295
1296
1297
1298
1299 /*
1300
1301
1302 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1303 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1304 AliXRDPROOFtoolkit tool;
1305 TChain * chainTr = tool.MakeChain("align.txt","Track",0,10200);
1306 chainTr->Lookup();
1307
1308 TCut c1pt("abs((tp1.fP[4]+tp2.fP[4])*0.5)<3"); // pt cut  - OK
1309 TCut cdy("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1310 TCut cdz("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1311 TCut cdphi("abs(tp1.fP[2]-tp2.fP[2])<0.02");
1312 TCut cdt("abs(tp1.fP[3]-tp2.fP[3])<0.02");
1313 TCut cd1pt("abs(tp1.fP[4]-tp2.fP[4])<0.3");    // delta 1/pt cut  -OK
1314 //
1315 //
1316 TCut cutA =  c1pt+cdy+cdz+cdphi+cdt+cd1pt;
1317 chainTr->Draw(">>listEL",acut,"entryList");
1318 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1319 chainTr->SetEntryList(elist);
1320
1321
1322
1323
1324
1325
1326 TCut cutS("s1%36==s2%36");
1327
1328 TCut cutN("c1>32&&c2>60");
1329 TCut cutC0("sqrt(tp2.fC[0])<1");
1330
1331 TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.4");
1332 TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.01");
1333 TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.01");
1334 TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
1335 TCut cutP=cutP0+cutP2+cutP3+cutP4+cutC0;
1336
1337 TCut cutX("abs(tp2.fX-133.6)<2");
1338
1339 TCut cutA = cutP+cutN;
1340
1341
1342 TCut cutY("abs(vcY.fElements-vtY.fElements)<0.3&&vcY.fElements!=0")
1343 TCut cutZ("abs(vcZ.fElements-vtZ.fElements)<0.3&&vcZ.fElements!=0")
1344
1345 */