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