1. Adding the AliExternalComparison
[u/mrichter/AliRoot.git] / 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 "TTree.h"
116 #include "TF1.h"
117 #include "TGraphErrors.h"
118 #include "AliTPCclusterMI.h"
119 #include "AliTPCseed.h"
120 #include "AliTracker.h"
121 #include "TClonesArray.h"
122 #include "AliExternalComparison.h"
123
124
125 #include "TTreeStream.h"
126 #include <iostream>
127 #include <sstream>
128 using namespace std;
129
130 AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0;
131 ClassImp(AliTPCcalibAlign)
132
133
134
135
136 AliTPCcalibAlign* AliTPCcalibAlign::Instance()
137 {
138   //
139   // Singleton implementation
140   // Returns an instance of this class, it is created if neccessary
141   //
142   if (fgInstance == 0){
143     fgInstance = new AliTPCcalibAlign();
144   }
145   return fgInstance;
146 }
147
148
149
150
151 AliTPCcalibAlign::AliTPCcalibAlign()
152   :  AliTPCcalibBase(),
153      fDphiHistArray(72*72),
154      fDthetaHistArray(72*72),
155      fDyHistArray(72*72),
156      fDzHistArray(72*72),
157      //
158      fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
159      fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
160      fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
161      fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
162      fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
163      fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     
164      fFitterArray12(72*72),
165      fFitterArray9(72*72),
166      fFitterArray6(72*72),
167      fMatrixArray12(72*72),
168      fMatrixArray9(72*72),
169      fMatrixArray6(72*72),
170      fCombinedMatrixArray6(72),
171      fCompTracklet(0),             // tracklet comparison
172      fNoField(kFALSE)
173 {
174   //
175   // Constructor
176   //
177   for (Int_t i=0;i<72*72;++i) {
178     fPoints[i]=0;
179   }
180 }
181
182 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
183   :AliTPCcalibBase(),  
184    fDphiHistArray(72*72),
185    fDthetaHistArray(72*72),
186    fDyHistArray(72*72),
187    fDzHistArray(72*72),
188    fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
189    fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
190    fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
191    fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
192    fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
193    fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     //
194    fFitterArray12(72*72),
195    fFitterArray9(72*72),
196    fFitterArray6(72*72),
197    fMatrixArray12(72*72),
198    fMatrixArray9(72*72),
199    fMatrixArray6(72*72),
200    fCombinedMatrixArray6(72),
201    fCompTracklet(0),             // tracklet comparison
202    fNoField(kFALSE)
203 {
204   //
205   // Constructor
206   //
207   SetName(name);
208   SetTitle(title);
209   for (Int_t i=0;i<72*72;++i) {
210     fPoints[i]=0;
211   }
212 }
213
214
215 AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
216   :AliTPCcalibBase(align),  
217    fDphiHistArray(align.fDphiHistArray),
218    fDthetaHistArray(align.fDthetaHistArray),
219    fDyHistArray(align.fDyHistArray),
220    fDzHistArray(align.fDzHistArray),
221    fDyPhiHistArray(align.fDyPhiHistArray),      // array of residual histograms  y     -kYPhi
222    fDzThetaHistArray(align.fDzThetaHistArray),    // array of residual histograms  z-z   -kZTheta
223    fDphiZHistArray(align.fDphiZHistArray),      // array of residual histograms  phi   -kPhiz
224    fDthetaZHistArray(align.fDthetaZHistArray),    // array of residual histograms  theta -kThetaz
225    fDyZHistArray(align.fDyZHistArray),        // array of residual histograms  y     -kYz
226    fDzZHistArray(align.fDzZHistArray),        // array of residual histograms  z     -kZz     
227    //
228    fFitterArray12(align.fFitterArray12),
229    fFitterArray9(align.fFitterArray9),
230    fFitterArray6(align.fFitterArray6),
231    
232    fMatrixArray12(align.fMatrixArray12),
233    fMatrixArray9(align.fMatrixArray9),
234    fMatrixArray6(align.fMatrixArray6),
235    fCombinedMatrixArray6(align.fCombinedMatrixArray6),
236    fCompTracklet(align.fCompTracklet),             // tracklet comparison
237    fNoField(align.fNoField)   
238 {
239   //
240   // copy constructor - copy also the content
241   //
242   TH1 * his = 0;
243   TObjArray * arr0=0;
244   const TObjArray *arr1=0;
245   for (Int_t index =0; index<72*72; index++){
246     for (Int_t iarray=0;iarray<10; iarray++){
247       if (iarray==kY){
248         arr0 = &fDyHistArray;
249         arr1 = &align.fDyHistArray;
250       }
251       if (iarray==kZ){
252         arr0 = &fDzHistArray;
253         arr1 = &align.fDzHistArray;
254       }
255       if (iarray==kPhi){
256         arr0 = &fDphiHistArray;
257         arr1 = &align.fDphiHistArray;
258       }
259       if (iarray==kTheta){
260         arr0 = &fDthetaHistArray;
261         arr1 = &align.fDthetaHistArray;
262       }
263       if (iarray==kYz){
264         arr0 = &fDyZHistArray;
265         arr1 = &align.fDyZHistArray;
266       }
267       if (iarray==kZz){
268         arr0 = &fDzZHistArray;
269         arr1 = &align.fDzZHistArray;
270       }
271       if (iarray==kPhiZ){
272         arr0 = &fDphiZHistArray;
273         arr1 = &align.fDphiZHistArray;
274       }
275       if (iarray==kThetaZ){
276         arr0 = &fDthetaZHistArray;
277         arr1 = &align.fDthetaZHistArray;
278       }
279
280       if (iarray==kYPhi){
281         arr0 = &fDyPhiHistArray;
282         arr1 = &align.fDyPhiHistArray;
283       }
284       if (iarray==kZTheta){
285         arr0 = &fDzThetaHistArray;
286         arr1 = &align.fDzThetaHistArray;
287       }
288
289       if (arr1->At(index)) {
290         his = (TH1*)arr1->At(index)->Clone();
291         his->SetDirectory(0);
292         arr0->AddAt(his,index);
293       }    
294     }
295   }
296 }
297
298
299 AliTPCcalibAlign::~AliTPCcalibAlign() {
300   //
301   // destructor
302   //
303 }
304
305 void AliTPCcalibAlign::Process(AliTPCseed *seed) {
306   //
307   // 
308   //
309   
310   TObjArray tracklets=
311     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
312                                     kFALSE,20,2);
313   //  TObjArray trackletsL=
314 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
315 //                                  kFALSE,20,2);
316 //   TObjArray trackletsQ=
317 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
318 //                                  kFALSE,20,2);
319 //   TObjArray trackletsR=
320 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
321 //                                  kFALSE,20,2);
322   tracklets.SetOwner();
323  //  trackletsL.SetOwner();
324 //   trackletsQ.SetOwner();
325 //   trackletsR.SetOwner();
326   if (tracklets.GetEntries()==2) {
327     AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
328     AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
329     if (t1->GetSector()>t2->GetSector()) {
330       AliTPCTracklet* tmp=t1;
331       t1=t2;
332       t2=tmp;
333     }
334     AliExternalTrackParam *common1=0,*common2=0;
335     if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
336       ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
337     delete common1;
338     delete common2;
339   }
340   
341 }
342
343 void AliTPCcalibAlign::Analyze(){
344   //
345   // Analyze function 
346   //
347   EvalFitters();
348 }
349
350
351 void AliTPCcalibAlign::Terminate(){
352   //
353   // Terminate function
354   // call base terminate + Eval of fitters
355   //
356   Info("AliTPCcalibAlign","Terminate");
357   EvalFitters();
358   AliTPCcalibBase::Terminate();
359 }
360
361
362
363
364 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
365                                         const AliExternalTrackParam &tp2,
366                                         const AliTPCseed * seed,
367                                         Int_t s1,Int_t s2) {
368
369   //
370   //
371   //
372   //
373   //
374   // Process function to fill fitters
375   //
376   Double_t t1[10],t2[10];
377   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
378   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
379   x1   =tp1.GetX();
380   y1   =tp1.GetY();
381   z1   =tp1.GetZ();
382   Double_t snp1=tp1.GetSnp();
383   dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
384   Double_t tgl1=tp1.GetTgl();
385   // dz/dx = 1/(cos(theta)*cos(phi))
386   dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1);
387   x2   =tp2.GetX();
388   y2   =tp2.GetY();
389   z2   =tp2.GetZ();
390   Double_t snp2=tp2.GetSnp();
391   dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
392   Double_t tgl2=tp2.GetTgl();
393   dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2);
394   Int_t accept =   AcceptTracklet(tp1,tp2);  
395   //
396   //
397   //
398   if (fStreamLevel>1 && seed){
399     TTreeSRedirector *cstream = GetDebugStreamer();
400     if (cstream){
401       static TVectorD vec1(5);
402       static TVectorD vec2(5);
403       vec1.SetElements(t1);
404       vec2.SetElements(t2);
405       AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
406       AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
407       (*cstream)<<"Tracklet"<<
408         "accept="<<accept<<
409         "run="<<fRun<<              //  run number
410         "event="<<fEvent<<          //  event number
411         "time="<<fTime<<            //  time stamp of event
412         "trigger="<<fTrigger<<      //  trigger
413         "mag="<<fMagF<<             //  magnetic field
414         "isOK="<<accept<<           //  flag - used for alignment
415         "tp1.="<<p1<<
416         "tp2.="<<p2<<
417         "v1.="<<&vec1<<
418         "v2.="<<&vec2<<
419         "s1="<<s1<<
420         "s2="<<s2<<
421         "\n";
422     }
423   }
424   if (accept>0) return;
425   t1[0]-=134.;
426   t2[0]-=134.;
427   t1[5]=0; t2[5]=0;
428   t1[6]=TMath::Sqrt(tp1.GetSigmaY2()+tp2.GetSigmaY2()); t2[6]=t1[6];
429   t1[7]=TMath::Sqrt(tp1.GetSigmaZ2()+tp2.GetSigmaZ2()); t2[7]=t1[7];
430   t1[8]=TMath::Sqrt(tp1.GetSigmaSnp2()+tp2.GetSigmaSnp2()); t2[8]=t1[8];
431   t1[9]=TMath::Sqrt(tp1.GetSigmaTgl2()+tp2.GetSigmaTgl2()); t2[9]=t1[9];
432
433   if (GetDebugLevel()>50) printf("Process track\n");
434   if (GetDebugLevel()>50) printf("Filling track\n");
435   //
436   // fill resolution histograms - previous cut included
437   if (seed) ProcessDiff(tp1,tp2, seed,s1,s2);
438   FillHisto(tp1,tp2,s1,s2);  
439   ProcessAlign(t1,t2,s1,s2);
440 }
441
442 void AliTPCcalibAlign::ProcessAlign(Double_t * t1,
443                                     Double_t * t2,
444                                     Int_t s1,Int_t s2){
445   //
446   // Do intersector alignment
447   //
448   Process12(t1,t2,GetOrMakeFitter12(s1,s2));
449   Process9(t1,t2,GetOrMakeFitter9(s1,s2));
450   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
451   ++fPoints[GetIndex(s1,s2)];
452 }
453
454 void AliTPCcalibAlign::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){
455   //
456   // Process the debug streamer tree
457   // Possible to modify selection criteria
458   // Used with entry list
459   //
460   TTreeSRedirector * cstream = new TTreeSRedirector("aligndump.root");
461
462   AliTPCcalibAlign *align = this;
463   //
464   TVectorD * vec1 = 0;
465   TVectorD * vec2 = 0;
466   AliExternalTrackParam * tp1 = 0;
467   AliExternalTrackParam * tp2 = 0;  
468   Int_t      s1 = 0;
469   Int_t      s2 = 0;                            
470   Int_t npoints =0;
471   {
472     Int_t entries=chainTracklet->GetEntries();
473     for (Int_t i=0; i< entries; i++){
474       chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
475       chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
476       chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
477       chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
478       chainTracklet->GetBranch("s1")->SetAddress(&s1);
479       chainTracklet->GetBranch("s2")->SetAddress(&s2);      
480       chainTracklet->GetEntry(i);
481       if (!vec1) continue;
482       if (!vec2) continue;
483       if (!tp1) continue;
484       if (!tp2) continue;
485       if (!vec1->GetMatrixArray()) continue;
486       if (!vec2->GetMatrixArray()) continue;
487       // make a local copy
488       AliExternalTrackParam par1(*tp1);
489       AliExternalTrackParam par2(*tp2);
490       TVectorD svec1(*vec1);
491       TVectorD svec2(*vec2);
492       //
493       if (s1==s2) continue;
494       if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i, npoints,s1,s2);
495       AliExternalTrackParam  cpar1(par1);
496       AliExternalTrackParam  cpar2(par2);      
497       Constrain1Pt(cpar1,par2,fNoField);
498       Constrain1Pt(cpar2,par1,fNoField);
499       Bool_t acceptComp = kFALSE;
500       if (comp) acceptComp=comp->AcceptPair(&par1,&par2);
501       if (comp) acceptComp&=comp->AcceptPair(&cpar1,&cpar2);
502       //
503       Int_t reject =   align->AcceptTracklet(par1,par2);
504       Int_t rejectC =align->AcceptTracklet(cpar1,cpar2); 
505
506       if (1||fStreamLevel>0){
507         (*cstream)<<"Tracklet"<<
508           "s1="<<s1<<
509           "s2="<<s2<<
510           "reject="<<reject<<
511           "rejectC="<<rejectC<<
512           "acceptComp="<<acceptComp<<
513           "tp1.="<<&par1<<
514           "tp2.="<<&par2<<      
515           "ctp1.="<<&cpar1<<
516           "ctp2.="<<&cpar2<<
517           "v1.="<<&svec1<<
518           "v2.="<<&svec2<<
519           "\n";
520       }
521       //
522       if (fNoField){
523         //
524         //
525       }
526       if (acceptComp) comp->Process(&cpar1,&cpar2);
527       //
528       if (reject>0 || rejectC>0) continue;
529       npoints++;
530       align->ProcessTracklets(cpar1,cpar2,0,s1,s2);
531       align->ProcessTracklets(cpar2,cpar1,0,s2,s1); 
532     }
533   }
534   delete cstream;
535 }
536
537
538 Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
539                                        const AliExternalTrackParam &p2){
540
541   //
542   // Accept pair of tracklets?
543   //
544   /*
545   // resolution cuts
546   TCut cutS0("sqrt(tp2.fC[0]+tp1.fC[0])<0.2");
547   TCut cutS1("sqrt(tp2.fC[2]+tp1.fC[2])<0.2");
548   TCut cutS2("sqrt(tp2.fC[5]+tp1.fC[5])<0.01");
549   TCut cutS3("sqrt(tp2.fC[9]+tp1.fC[9])<0.01");
550   TCut cutS4("sqrt(tp2.fC[14]+tp1.fC[14])<0.25");
551   TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
552   //
553   // parameters matching cuts
554   TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
555   TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
556   TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03");
557   TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03");
558   TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
559   TCut cutPP4("abs(tp1.fP[4]-tp2.fP[4])/sqrt(tp2.fC[14]+tp1.fC[14])<3");
560   TCut cutP=cutP0+cutP1+cutP2+cutP3+cutP4+cutPP4;
561   */  
562   //
563   // resolution cuts
564   Int_t reject=0;
565   const Double_t *cp1 = p1.GetCovariance();
566   const Double_t *cp2 = p2.GetCovariance();
567   if (TMath::Sqrt(cp1[0]+cp2[0])>0.2)  reject|=1;;
568   if (TMath::Sqrt(cp1[2]+cp2[2])>0.2)  reject|=2;
569   if (TMath::Sqrt(cp1[5]+cp2[5])>0.01) reject|=4;
570   if (TMath::Sqrt(cp1[9]+cp2[9])>0.01) reject|=8;
571   if (TMath::Sqrt(cp1[14]+cp2[14])>0.2) reject|=16;
572
573   //parameters difference
574   const Double_t *tp1 = p1.GetParameter();
575   const Double_t *tp2 = p2.GetParameter();
576   if (TMath::Abs(tp1[0]-tp2[0])>0.6) reject|=32;
577   if (TMath::Abs(tp1[1]-tp2[1])>0.6) reject|=64;
578   if (TMath::Abs(tp1[2]-tp2[2])>0.03) reject|=128;
579   if (TMath::Abs(tp1[3]-tp2[3])>0.03) reject|=526;
580   if (TMath::Abs(tp1[4]-tp2[4])>0.4) reject|=1024;
581   if (TMath::Abs(tp1[4]-tp2[4])/TMath::Sqrt(cp1[14]+cp2[14])>4) reject|=2048;
582   
583   //
584   if (TMath::Abs(tp2[1])>235) reject|=2*4096;
585   
586   if (fNoField){
587     
588   }
589
590   return reject;
591 }
592
593
594
595 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
596                                     const AliExternalTrackParam &t2,
597                                     const AliTPCseed *seed,
598                                     Int_t s1,Int_t s2)
599 {
600   //
601   // Process local residuals function
602   // 
603   TVectorD vecX(160);
604   TVectorD vecY(160);
605   TVectorD vecZ(160);
606   TVectorD vecClY(160);
607   TVectorD vecClZ(160);
608   TClonesArray arrCl("AliTPCclusterMI",160);
609   arrCl.ExpandCreateFast(160);
610   Int_t count1=0, count2=0;
611   
612   for (Int_t i=0;i<160;++i) {
613     AliTPCclusterMI *c=seed->GetClusterPointer(i);
614     vecX[i]=0;
615     vecY[i]=0;
616     vecZ[i]=0;
617     if (!c) continue;
618     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
619     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
620     vecClY[i] = c->GetY();
621     vecClZ[i] = c->GetZ();
622     cl=*c;
623     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
624     if (c->GetDetector()==s1) ++count1;
625     if (c->GetDetector()==s2) ++count2;
626     Double_t gxyz[3],xyz[3];
627     t1.GetXYZ(gxyz);
628     Float_t bz = AliTracker::GetBz(gxyz);
629     par->GetYAt(c->GetX(), bz, xyz[1]);
630     par->GetZAt(c->GetX(), bz, xyz[2]);
631     vecX[i] = c->GetX();
632     vecY[i]= xyz[1];
633     vecZ[i]= xyz[2];
634   }
635   //
636   //
637   if (fStreamLevel>5){
638     //
639     // huge output - cluster residuals to be investigated
640     //
641     TTreeSRedirector *cstream = GetDebugStreamer();
642     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
643     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
644     /*
645       
646       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");
647
648     */
649
650     if (cstream){
651       (*cstream)<<"Track"<<
652         "run="<<fRun<<              //  run number
653         "event="<<fEvent<<          //  event number
654         "time="<<fTime<<            //  time stamp of event
655         "trigger="<<fTrigger<<      //  trigger
656         "mag="<<fMagF<<             //  magnetic field
657         "Cl.="<<&arrCl<<
658         //"tp0.="<<p0<<
659         "tp1.="<<p1<<
660         "tp2.="<<p2<<
661         "vtX.="<<&vecX<<
662         "vtY.="<<&vecY<<
663         "vtZ.="<<&vecZ<<
664         "vcY.="<<&vecClY<<
665         "vcZ.="<<&vecClZ<<
666         "s1="<<s1<<
667         "s2="<<s2<<
668         "c1="<<count1<<
669         "c2="<<count2<<
670         "\n";
671     }
672   }
673 }
674
675
676
677
678 void AliTPCcalibAlign::Process12(const Double_t *t1,
679                                  const Double_t *t2,
680                                  TLinearFitter *fitter) {
681   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
682   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
683   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
684   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
685   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
686   //
687   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
688   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
689   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
690
691
692
693   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
694   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
695
696   //
697   Double_t sy    = t1[6];
698   Double_t sz    = t1[7];
699   Double_t sdydx = t1[8];
700   Double_t sdzdx = t1[9];
701
702   Double_t p[12];
703   Double_t value;
704
705   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
706   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
707   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
708   for (Int_t i=0; i<12;i++) p[i]=0.;
709   p[3+0] = x1;          // a10
710   p[3+1] = y1;          // a11
711   p[3+2] = z1;          // a12
712   p[9+1] = 1.;          // a13
713   p[0+1] = y1*dydx2;    // a01
714   p[0+2] = z1*dydx2;    // a02
715   p[9+0] = dydx2;       // a03
716   value  = y2;
717   fitter->AddPoint(p,value,sy);
718
719   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
720   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
721   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
722   for (Int_t i=0; i<12;i++) p[i]=0.;
723   p[6+0] = x1;           // a20 
724   p[6+1] = y1;           // a21
725   p[6+2] = z1;           // a22
726   p[9+2] = 1.;           // a23
727   p[0+1] = y1*dzdx2;     // a01
728   p[0+2] = z1*dzdx2;     // a02
729   p[9+0] = dzdx2;        // a03
730   value  = z2;
731   fitter->AddPoint(p,value,sz);
732
733   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
734   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
735   for (Int_t i=0; i<12;i++) p[i]=0.;
736   p[3+0] = 1.;           // a10
737   p[3+1] = dydx1;        // a11
738   p[3+2] = dzdx1;        // a12
739   p[0+0] = -dydx2;       // a00
740   p[0+1] = -dydx1*dydx2; // a01
741   p[0+2] = -dzdx1*dydx2; // a02
742   value  = 0.;
743   fitter->AddPoint(p,value,sdydx);
744
745   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
746   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
747   for (Int_t i=0; i<12;i++) p[i]=0.;
748   p[6+0] = 1;            // a20
749   p[6+1] = dydx1;        // a21
750   p[6+2] = dzdx1;        // a22
751   p[0+0] = -dzdx2;       // a00
752   p[0+1] = -dydx1*dzdx2; // a01
753   p[0+2] = -dzdx1*dzdx2; // a02
754   value  = 0.;
755   fitter->AddPoint(p,value,sdzdx);
756 }
757
758 void AliTPCcalibAlign::Process9(Double_t *t1,
759                                 Double_t *t2,
760                                 TLinearFitter *fitter) {
761   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
762   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
763   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
764   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
765   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
766   //
767   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
768   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
769   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8] 
770
771
772   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
773   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
774   //
775   Double_t sy    = t1[6];
776   Double_t sz    = t1[7];
777   Double_t sdydx = t1[8];
778   Double_t sdzdx = t1[9];
779   //
780   Double_t p[12];
781   Double_t value;
782
783   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
784   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
785   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
786   for (Int_t i=0; i<12;i++) p[i]=0.;
787   p[2]   += x1;           // a10
788   //p[]  +=1;             // a11
789   p[3]   += z1;           // a12    
790   p[7]   += 1;            // a13
791   p[0]   += y1*dydx2;     // a01
792   p[1]   += z1*dydx2;     // a02
793   p[6]   += dydx2;        // a03
794   value   = y2-y1;        //-a11
795   fitter->AddPoint(p,value,sy);
796   //
797   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
798   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
799   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
800   for (Int_t i=0; i<12;i++) p[i]=0.;
801   p[4]   += x1;           // a20 
802   p[5]   += y1;           // a21
803   //p[]  += z1;           // a22
804   p[8]   += 1.;           // a23
805   p[0]   += y1*dzdx2;     // a01
806   p[1]   += z1*dzdx2;     // a02
807   p[6]   += dzdx2;        // a03
808   value  = z2-z1;         //-a22
809   fitter->AddPoint(p,value,sz);
810
811   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
812   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
813   for (Int_t i=0; i<12;i++) p[i]=0.;
814   p[2]   += 1.;           // a10
815   //p[]  += dydx1;      // a11
816   p[3]   += dzdx1;        // a12
817   //p[]  += -dydx2;       // a00
818   p[0]   += -dydx1*dydx2; // a01
819   p[1]   += -dzdx1*dydx2; // a02
820   value  = -dydx1+dydx2;  // -a11 + a00
821   fitter->AddPoint(p,value,sdydx);
822
823   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
824   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
825   for (Int_t i=0; i<12;i++) p[i]=0.;
826   p[4]   += 1;            // a20
827   p[5]   += dydx1;        // a21
828   //p[]  += dzdx1;        // a22
829   //p[]  += -dzdx2;       // a00
830   p[0]   += -dydx1*dzdx2; // a01
831   p[1]   += -dzdx1*dzdx2; // a02
832   value  = -dzdx1+dzdx2;  // -a22 + a00
833   fitter->AddPoint(p,value,sdzdx);
834 }
835
836 void AliTPCcalibAlign::Process6(Double_t *t1,
837                                 Double_t *t2,
838                                 TLinearFitter *fitter) {
839   // x2    =  1  *x1 +-a01*y1 + 0      +a03
840   // y2    =  a01*x1 + 1  *y1 + 0      +a13
841   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
842   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
843   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
844   //
845   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
846   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
847   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5] 
848
849   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
850   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
851
852   //
853   Double_t sy    = t1[6];
854   Double_t sz    = t1[7];
855   Double_t sdydx = t1[8];
856   Double_t sdzdx = t1[9];
857
858   Double_t p[12];
859   Double_t value;
860   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
861   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
862   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
863   for (Int_t i=0; i<12;i++) p[i]=0.;
864   p[0]   += x1;           // a10
865   //p[]  +=1;             // a11
866   //p[]  += z1;           // a12    
867   p[4]   += 1;            // a13
868   p[0]   += -y1*dydx2;    // a01
869   //p[]  += z1*dydx2;     // a02
870   p[3]   += dydx2;        // a03
871   value   = y2-y1;        //-a11
872   fitter->AddPoint(p,value,sy);
873   //
874   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
875   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
876   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
877   for (Int_t i=0; i<12;i++) p[i]=0.;
878   p[1]   += x1;           // a20 
879   p[2]   += y1;           // a21
880   //p[]  += z1;           // a22
881   p[5]   += 1.;           // a23
882   p[0]   += -y1*dzdx2;    // a01
883   //p[]   += z1*dzdx2;     // a02
884   p[3]   += dzdx2;        // a03
885   value  = z2-z1;         //-a22
886   fitter->AddPoint(p,value,sz);
887
888   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
889   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
890   for (Int_t i=0; i<12;i++) p[i]=0.;
891   p[0]   += 1.;           // a10
892   //p[]  += dydx1;      // a11       
893   //p[]   += dzdx1;        // a12
894   //p[]  += -dydx2;       // a00
895   //p[0]   +=  dydx1*dydx2; // a01         FIXME- 0912 MI
896   //p[]   += -dzdx1*dydx2; // a02
897   value  = -dydx1+dydx2;  // -a11 + a00
898   fitter->AddPoint(p,value,sdydx);
899
900   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
901   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
902   for (Int_t i=0; i<12;i++) p[i]=0.;
903   p[1]   += 1;            // a20
904   //  p[2]   += dydx1;        // a21   FIXME- 0912 MI
905   //p[]  += dzdx1;        // a22
906   //p[]  += -dzdx2;       // a00
907   //p[0]   +=  dydx1*dzdx2; // a01     FIXME- 0912 MI
908   //p[]  += -dzdx1*dzdx2; // a02
909   value  = -dzdx1+dzdx2;  // -a22 + a00
910   fitter->AddPoint(p,value,sdzdx);
911 }
912
913
914
915
916 void AliTPCcalibAlign::EvalFitters() {
917   //
918   // Analyze function 
919   // 
920   // Perform the fitting using linear fitters
921   //
922   Int_t kMinPoints =50;
923   TLinearFitter *f;
924   TFile fff("alignDebug.root","recreate");
925   for (Int_t s1=0;s1<72;++s1)
926     for (Int_t s2=0;s2<72;++s2){
927       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
928         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
929         if (f->Eval()!=0) {
930           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
931           f->Write(Form("f12_%d_%d",s1,s2));
932         }else{
933           f->Write(Form("f12_%d_%d",s1,s2));
934         }
935       }
936       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
937         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
938         if (f->Eval()!=0) {
939           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
940         }else{
941           f->Write(Form("f9_%d_%d",s1,s2));
942         }
943       }
944       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
945         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
946         if (f->Eval()!=0) {
947           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
948         }else{
949           f->Write(Form("f6_%d_%d",s1,s2));
950         }
951       }
952     }
953   TMatrixD mat(4,4);
954   for (Int_t s1=0;s1<72;++s1)
955     for (Int_t s2=0;s2<72;++s2){
956       if (GetTransformation12(s1,s2,mat)){
957         fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2));
958       }
959       if (GetTransformation9(s1,s2,mat)){
960         fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2));
961       }
962       if (GetTransformation6(s1,s2,mat)){
963         fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2));
964       }
965     }
966   //this->Write("align");
967   
968 }
969
970 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
971   //
972   // get or make fitter - general linear transformation
973   //
974   static Int_t counter12=0;
975   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]");
976   TLinearFitter * fitter = GetFitter12(s1,s2);
977   if (fitter) return fitter;
978   //  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]");
979   fitter =new TLinearFitter(&f12,"");
980   fitter->StoreData(kFALSE);
981   fFitterArray12.AddAt(fitter,GetIndex(s1,s2)); 
982   counter12++;
983   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
984   return fitter;
985 }
986
987 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
988   //
989   //get or make fitter - general linear transformation - no scaling
990   // 
991   static Int_t counter9=0;
992   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
993   TLinearFitter * fitter = GetFitter9(s1,s2);
994   if (fitter) return fitter;
995   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
996   fitter =new TLinearFitter(&f9,"");
997   fitter->StoreData(kFALSE);
998   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
999   counter9++;
1000   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
1001   return fitter;
1002 }
1003
1004 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
1005   //
1006   // get or make fitter  - 6 paramater linear tranformation
1007   //                     - no scaling
1008   //                     - rotation x-y
1009   //                     - tilting x-z, y-z
1010   static Int_t counter6=0;
1011   static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1012   TLinearFitter * fitter = GetFitter6(s1,s2);
1013   if (fitter) return fitter;
1014   //  fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1015   fitter=new TLinearFitter(&f6,"");
1016   fitter->StoreData(kTRUE);
1017   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
1018   counter6++;
1019   if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
1020   return fitter;
1021 }
1022
1023
1024
1025
1026
1027 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
1028   //
1029   // GetTransformation matrix - 12 paramaters - generael linear transformation
1030   //
1031   if (!GetFitter12(s1,s2))
1032     return false;
1033   else {
1034     TVectorD p(12);
1035     GetFitter12(s1,s2)->GetParameters(p);
1036     a.ResizeTo(4,4);
1037     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
1038     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
1039     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
1040     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1041     return true;
1042   } 
1043 }
1044
1045 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
1046   //
1047   // GetTransformation matrix - 9 paramaters - general linear transformation
1048   //                            No scaling
1049   //
1050   if (!GetFitter9(s1,s2))
1051     return false;
1052   else {
1053     TVectorD p(9);
1054     GetFitter9(s1,s2)->GetParameters(p);
1055     a.ResizeTo(4,4);
1056     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
1057     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
1058     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
1059     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1060     return true;
1061   } 
1062 }
1063
1064 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
1065   //
1066   // GetTransformation matrix - 6  paramaters
1067   //                            3  translation
1068   //                            1  rotation -x-y  
1069   //                            2  tilting x-z y-z
1070   if (!GetFitter6(s1,s2))
1071     return false;
1072   else {
1073     TVectorD p(6);
1074     GetFitter6(s1,s2)->GetParameters(p);
1075     a.ResizeTo(4,4);
1076     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
1077     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
1078     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
1079     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1080     return true;
1081   } 
1082 }
1083
1084 void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
1085                                         const AliExternalTrackParam &tp2,
1086                                         Int_t s1,Int_t s2) {
1087   //
1088   // Fill residual histograms
1089   // Innner-Outer
1090   // Left right - x-y
1091   // A-C side
1092   if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0)  {  
1093     GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));    
1094     GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1095     GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
1096     GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
1097     //
1098     GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));    
1099     GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1100     GetHisto(kYz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetY()-tp2.GetY());
1101     GetHisto(kZz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetZ()-tp2.GetZ());
1102     //
1103     GetHisto(kYPhi,s1,s2,kTRUE)->Fill(tp1.GetSnp(),tp1.GetY()-tp2.GetY());
1104     GetHisto(kZTheta,s1,s2,kTRUE)->Fill(tp1.GetTgl(),tp1.GetZ()-tp2.GetZ());
1105
1106
1107   }  
1108 }
1109
1110
1111
1112 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
1113 {
1114   //
1115   // return specified residual histogram - it is only QA
1116   // if force specified the histogram and given histogram is not existing 
1117   //  new histogram is created
1118   //
1119   if (GetIndex(s1,s2)>=72*72) return 0;
1120   TObjArray *histoArray=0;
1121   switch (type) {
1122   case kY:
1123     histoArray = &fDyHistArray; break;
1124   case kZ:
1125     histoArray = &fDzHistArray; break;
1126   case kPhi:
1127     histoArray = &fDphiHistArray; break;
1128   case kTheta:
1129     histoArray = &fDthetaHistArray; break;
1130   case kYPhi:
1131     histoArray = &fDyPhiHistArray; break;
1132   case kZTheta:
1133     histoArray = &fDzThetaHistArray; break;
1134   case kYz:
1135     histoArray = &fDyZHistArray; break;
1136   case kZz:
1137     histoArray = &fDzZHistArray; break;
1138   case kPhiZ:
1139     histoArray = &fDphiZHistArray; break;
1140   case kThetaZ:
1141     histoArray = &fDthetaZHistArray; break;
1142   }
1143   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
1144   if (histo) return histo;
1145   if (force==kFALSE) return 0; 
1146   //
1147   stringstream name;
1148   stringstream title;
1149   switch (type) {    
1150   case kY:
1151     name<<"hist_y_"<<s1<<"_"<<s2;
1152     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1153     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1154     break;
1155   case kZ:
1156     name<<"hist_z_"<<s1<<"_"<<s2;
1157     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1158     histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1159     break;
1160   case kPhi:
1161     name<<"hist_phi_"<<s1<<"_"<<s2;
1162     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1163     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1164     break;
1165   case kTheta:
1166     name<<"hist_theta_"<<s1<<"_"<<s2;
1167     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1168     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1169     break;
1170     //
1171     //
1172   case kYPhi:
1173     name<<"hist_yphi_"<<s1<<"_"<<s2;
1174     title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1175     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,128,-0.3,0.3); // +/- 3 mm
1176     break;
1177   case kZTheta:
1178     name<<"hist_ztheta_"<<s1<<"_"<<s2;
1179     title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1180     histo = new TH2F(name.str().c_str(),title.str().c_str(),128,20,-1,1,-0.3,0.3); // +/- 3 mm
1181     break;
1182     //
1183     //
1184     //
1185   case kYz:
1186     name<<"hist_yz_"<<s1<<"_"<<s2;
1187     title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1188     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1189     break;
1190   case kZz:
1191     name<<"hist_zz_"<<s1<<"_"<<s2;
1192     title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1193     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1194     break;
1195   case kPhiZ:
1196     name<<"hist_phiz_"<<s1<<"_"<<s2;
1197     title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1198     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1199     break;
1200   case kThetaZ:
1201     name<<"hist_thetaz_"<<s1<<"_"<<s2;
1202     title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1203     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1204     break;
1205
1206
1207   }
1208   histo->SetDirectory(0);
1209   histoArray->AddAt(histo,GetIndex(s1,s2));
1210   return histo;
1211 }
1212
1213 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec, 
1214                                            Int_t i0, Int_t i1, FitType type) 
1215 {
1216   //
1217   //
1218   //
1219   TMatrixD mat;
1220   //TObjArray *fitArray=0;
1221   Double_t xsec[1000];
1222   Double_t ysec[1000];
1223   Int_t npoints=0;
1224   for (Int_t isec = sec0; isec<=sec1; isec++){
1225     Int_t isec2 = (isec+dsec)%72;    
1226     switch (type) {
1227     case k6:
1228       GetTransformation6(isec,isec2,mat);break;
1229     case k9:
1230       GetTransformation9(isec,isec2,mat);break;
1231     case k12:
1232       GetTransformation12(isec,isec2,mat);break;
1233     }
1234     xsec[npoints]=isec;
1235     ysec[npoints]=mat(i0,i1);
1236     ++npoints;
1237   }
1238   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1239   Char_t name[1000];
1240   sprintf(name,"Mat[%d,%d]  Type=%d",i0,i1,type);
1241   gr->SetName(name);
1242   return gr;
1243 }
1244
1245 void  AliTPCcalibAlign::MakeTree(const char *fname){
1246   //
1247   // make tree with alignment cosntant  -
1248   // For  QA visualization
1249   //
1250   /*
1251    TFile f("CalibObjects.root");
1252    TObjArray *array  = (TObjArray*)f.Get("TPCCalib");
1253    AliTPCcalibAlign   *alignTPC = (AliTPCcalibAlign   *)array->At(0);
1254    alignTPC->MakeTree("alignTree.root");
1255    TFile falign("alignTree.root");
1256    Align->Draw("dy")
1257    */
1258   const Int_t kMinPoints=50;
1259   TTreeSRedirector cstream(fname);
1260   for (Int_t s1=0;s1<72;++s1)
1261     for (Int_t s2=0;s2<72;++s2){
1262       TMatrixD m6;
1263       TMatrixD m6FX;
1264       TMatrixD m9;
1265       TMatrixD m12;
1266       Double_t dy=0, dz=0, dphi=0,dtheta=0;
1267       Double_t sy=0, sz=0, sphi=0,stheta=0;
1268       Double_t ny=0, nz=0, nphi=0,ntheta=0;
1269       Double_t chi2v12=0, chi2v9=0, chi2v6=0;
1270       Int_t npoints=0;
1271       TLinearFitter * fitter = 0;      
1272       if (fPoints[GetIndex(s1,s2)]>kMinPoints){
1273         //
1274         //
1275         //
1276         fitter = GetFitter12(s1,s2);
1277         npoints = fitter->GetNpoints();
1278         chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1279         //
1280         fitter = GetFitter9(s1,s2);
1281         npoints = fitter->GetNpoints();
1282         chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1283         //
1284         fitter = GetFitter6(s1,s2);
1285         npoints = fitter->GetNpoints();
1286         chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1287
1288         //
1289         GetTransformation6(s1,s2,m6);
1290         GetTransformation9(s1,s2,m9);
1291         GetTransformation12(s1,s2,m12);
1292         //
1293         fitter = GetFitter6(s1,s2);
1294         fitter->FixParameter(3,0);
1295         fitter->Eval();
1296         GetTransformation6(s1,s2,m6FX);
1297         //
1298         TH1 * his=0;
1299         his = GetHisto(kY,s1,s2);
1300         if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1301         his = GetHisto(kZ,s1,s2);
1302         if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1303         his = GetHisto(kPhi,s1,s2);
1304         if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
1305         his = GetHisto(kTheta,s1,s2);
1306         if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
1307         //
1308
1309       }
1310
1311       // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1312       // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1313       // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1314       // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1315       // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1316       //
1317       //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
1318       //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
1319       //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
1320       
1321       //
1322       // 
1323       // dy:-(134*m6.fElements[4]+m6.fElements[7])
1324       // 
1325       // dphi:-(m6.fElements[4])
1326       //
1327       // dz:134*m6.fElements[8]+m6.fElements[11]
1328       //
1329       // dtheta:m6.fElements[8]
1330       //
1331       cstream<<"Align"<<
1332         "s1="<<s1<<     // reference sector
1333         "s2="<<s2<<     // sector to align
1334         "m6FX.="<<&m6FX<<   // tranformation matrix
1335         "m6.="<<&m6<<   // tranformation matrix
1336         "m9.="<<&m9<<   // 
1337         "m12.="<<&m12<<
1338         "chi2v12="<<chi2v12<<
1339         "chi2v9="<<chi2v9<<
1340         "chi2v6="<<chi2v6<<
1341         //               histograms mean RMS and entries
1342         "dy="<<dy<<  
1343         "sy="<<sy<<
1344         "ny="<<ny<<
1345         "dz="<<dz<<
1346         "sz="<<sz<<
1347         "nz="<<nz<<
1348         "dphi="<<dphi<<
1349         "sphi="<<sphi<<
1350         "nphi="<<nphi<<
1351         "dtheta="<<dtheta<<
1352         "stheta="<<stheta<<
1353         "ntheta="<<ntheta<<
1354         "\n";
1355     }
1356
1357 }
1358
1359
1360 //_____________________________________________________________________
1361 Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
1362   //
1363   // merge function 
1364   //
1365   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
1366   if (!list)
1367     return 0;  
1368   if (list->IsEmpty())
1369     return 1;
1370   
1371   TIterator* iter = list->MakeIterator();
1372   TObject* obj = 0;  
1373   iter->Reset();
1374   Int_t count=0;
1375   //
1376   TString str1(GetName());
1377   while((obj = iter->Next()) != 0)
1378     {
1379       AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
1380       if (entry == 0) continue; 
1381       if (str1.CompareTo(entry->GetName())!=0) continue;
1382       Add(entry);
1383       count++;
1384     } 
1385   return count;
1386 }
1387
1388
1389 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
1390   //
1391   // Add entry - used for merging of compoents
1392   //
1393   for (Int_t i=0; i<72;i++){
1394     for (Int_t j=0; j<72;j++){
1395       if (align->fPoints[GetIndex(i,j)]<10) continue;
1396       fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
1397       //
1398       //
1399       //
1400       for (Int_t itype=0; itype<10; itype++){
1401         TH1 * his0=0, *his1=0;
1402         his0 = GetHisto((HistoType)itype,i,j);
1403         his1 = align->GetHisto((HistoType)itype,i,j);
1404         if (his1){
1405           if (his0) his0->Add(his1);
1406           else {
1407             his0 = GetHisto(kY,i,j,kTRUE);
1408             his0->Add(his1);
1409           }
1410         }       
1411       }      
1412     }
1413   }
1414   TLinearFitter *f0=0;
1415   TLinearFitter *f1=0;
1416   for (Int_t i=0; i<72;i++){
1417     for (Int_t j=0; j<72;j++){     
1418       if (align->fPoints[GetIndex(i,j)]<20) continue;
1419       //
1420       //
1421       // fitter12
1422       f0 =  GetFitter12(i,j);
1423       f1 =  align->GetFitter12(i,j);
1424       if (f1 &&f1->Eval()){
1425         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1426         else {
1427           fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
1428         }
1429       }      
1430       //
1431       // fitter9
1432       f0 =  GetFitter9(i,j);
1433       f1 =  align->GetFitter9(i,j);
1434       if (f1&&f1->Eval()){
1435         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1436         else { 
1437           fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
1438         }
1439       }      
1440       f0 =  GetFitter6(i,j);
1441       f1 =  align->GetFitter6(i,j);
1442       if (f1 &&f1->Eval()){
1443         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1444         else {
1445           fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
1446         }
1447       }   
1448     }
1449   }
1450 }
1451
1452 Double_t AliTPCcalibAlign::Correct(Int_t type, Int_t value, Int_t s1, Int_t s2, Double_t x1, Double_t y1, Double_t z1, Double_t dydx1,Double_t dzdx1){
1453   //
1454   // GetTransformed value
1455   //
1456   //
1457   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1458   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1459   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1460   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1461   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1462   
1463   
1464   const TMatrixD * mat = GetTransformation(s1,s2,type);
1465   if (!mat) {
1466     if (value==0) return x1;
1467     if (value==1) return y1;
1468     if (value==2) return z1;
1469     if (value==3) return dydx1;
1470     if (value==4) return dzdx1;
1471     //
1472     if (value==5) return dydx1;
1473     if (value==6) return dzdx1;
1474     return 0;
1475   }
1476   Double_t valT=0;
1477
1478   if (value==0){
1479     valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
1480   }
1481
1482   if (value==1){
1483     valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
1484   }
1485   if (value==2){
1486     valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
1487   }
1488   if (value==3){
1489     //    dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1490     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
1491     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1492   }
1493
1494   if (value==4){
1495     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)    
1496     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
1497     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1498   }
1499   //
1500   if (value==5){
1501     // onlys shift in angle
1502     //    dydx2 =  (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1503     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1;
1504   }
1505
1506   if (value==6){
1507     // only shift in angle
1508     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)    
1509     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1;
1510   }
1511   //
1512   return valT;
1513 }
1514
1515
1516 void  AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){
1517   //
1518   // Update track parameters t1
1519   //
1520   TMatrixD vecXk(5,1);    // X vector
1521   TMatrixD covXk(5,5);    // X covariance 
1522   TMatrixD matHk(1,5);    // vector to mesurement
1523   TMatrixD measR(1,1);    // measurement error 
1524   //TMatrixD matQk(5,5);    // prediction noise vector
1525   TMatrixD vecZk(1,1);    // measurement
1526   //
1527   TMatrixD vecYk(1,1);    // Innovation or measurement residual
1528   TMatrixD matHkT(5,1);
1529   TMatrixD matSk(1,1);    // Innovation (or residual) covariance
1530   TMatrixD matKk(5,1);    // Optimal Kalman gain
1531   TMatrixD mat1(5,5);     // update covariance matrix
1532   TMatrixD covXk2(5,5);   // 
1533   TMatrixD covOut(5,5);
1534   //
1535   Double_t *param1=(Double_t*) track1.GetParameter();
1536   Double_t *covar1=(Double_t*) track1.GetCovariance();
1537
1538   //
1539   // copy data to the matrix
1540   for (Int_t ipar=0; ipar<5; ipar++){
1541     vecXk(ipar,0) = param1[ipar];
1542     for (Int_t jpar=0; jpar<5; jpar++){
1543       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
1544     }
1545   }
1546   //
1547   //
1548   //
1549   vecZk(0,0) = track2.GetParameter()[4];   // 1/pt measurement from track 2
1550   measR(0,0) = track2.GetCovariance()[14];  // 1/pt measurement error
1551   if (noField) {
1552     measR(0,0)*=0.000000001;
1553     vecZk(0,0)=0.;
1554   }
1555   //
1556   matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0;  
1557   matHk(0,3)= 0;    matHk(0,4)= 1;           // vector to measurement
1558   //
1559   //
1560   //
1561   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
1562   matHkT=matHk.T(); matHk.T();
1563   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
1564   matSk.Invert();
1565   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
1566   vecXk += matKk*vecYk;                      //  updated vector 
1567   mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
1568   covXk2 = (mat1-(matKk*matHk));
1569   covOut =  covXk2*covXk; 
1570   //
1571   //
1572   //
1573   // copy from matrix to parameters
1574   if (0) {
1575     covOut.Print();
1576     vecXk.Print();
1577     covXk.Print();
1578     track1.Print();
1579     track2.Print();
1580   }
1581
1582   for (Int_t ipar=0; ipar<5; ipar++){
1583     param1[ipar]= vecXk(ipar,0) ;
1584     for (Int_t jpar=0; jpar<5; jpar++){
1585       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
1586     }
1587   }
1588
1589 }
1590
1591 void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){
1592   //
1593   //  Global Align -combine the partial alignment of pair of sectors
1594   //  minPoints - minimal number of points - don't use sector alignment wit smaller number
1595   //  sysError  - error added to the alignemnt error
1596   //
1597   AliTPCcalibAlign * align = this;
1598   TMatrixD * arrayAlign[72];
1599   TMatrixD * arrayAlignDiff[72];
1600   //
1601   for (Int_t i=0;i<72; i++) {
1602     TMatrixD * mat = new TMatrixD(4,4);
1603     mat->UnitMatrix();
1604     arrayAlign[i]=mat;
1605     arrayAlignDiff[i]=(TMatrixD*)(mat->Clone());
1606   }
1607
1608   TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root");
1609   for (Int_t iter=0; iter<niter;iter++){
1610     printf("Iter=\t%d\n",iter);
1611     for (Int_t is0=0;is0<72; is0++) {
1612       //
1613       //TMatrixD  *mati0 = arrayAlign[is0];
1614       TMatrixD matDiff(4,4);
1615       Double_t sumw=0;      
1616       for (Int_t is1=0;is1<72; is1++) {
1617         Bool_t invers=kFALSE;
1618         Int_t npoints=0;
1619         TMatrixD covar;
1620         TVectorD errors;
1621         const TMatrixD *mat = align->GetTransformation(is0,is1,0); 
1622         if (mat){
1623           npoints = align->GetFitter6(is0,is1)->GetNpoints();
1624           if (npoints>minPoints){
1625             align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar);
1626             align->GetFitter6(is0,is1)->GetErrors(errors);
1627           }
1628         }
1629         else{
1630           invers=kTRUE;
1631           mat = align->GetTransformation(is1,is0,0); 
1632           if (mat) {
1633             npoints = align->GetFitter6(is1,is0)->GetNpoints();
1634             if (npoints>minPoints){
1635               align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar);
1636               align->GetFitter6(is1,is0)->GetErrors(errors);
1637             }
1638           }
1639         }
1640         if (!mat) continue;
1641         if (npoints<minPoints) continue;
1642         //
1643         Double_t weight=1;
1644         if (is1/36>is0/36) weight*=2./3.; //IROC-OROC
1645         if (is1/36<is0/36) weight*=1./3.; //OROC-IROC
1646         if (is1/36==is0/36) weight*=1/3.; //OROC-OROC
1647         if (is1%36!=is0%36) weight*=1/2.; //Not up-down
1648         weight/=(errors[4]*errors[4]+sysError*sysError);
1649         //
1650         //
1651         TMatrixD matT = *mat;   
1652         if (invers)  matT.Invert();
1653         TMatrixD diffMat= (*(arrayAlign[is1]))*matT;
1654         diffMat-=(*arrayAlign[is0]);
1655         matDiff+=weight*diffMat;
1656         sumw+=weight;
1657
1658         (*cstream)<<"LAlign"<<
1659           "iter="<<iter<<
1660           "s0="<<is0<<
1661           "s1="<<is1<<
1662           "npoints="<<npoints<<
1663           "m60.="<<arrayAlign[is0]<<
1664           "m61.="<<arrayAlign[is1]<<
1665           "m01.="<<&matT<<
1666           "diff.="<<&diffMat<<
1667           "cov.="<<&covar<<
1668           "err.="<<&errors<<
1669           "\n";
1670       }
1671       if (sumw>0){
1672         matDiff*=1/sumw;
1673         matDiff(0,0)=0;
1674         matDiff(1,1)=0;
1675         matDiff(1,1)=0;
1676         matDiff(1,1)=0; 
1677         (*arrayAlignDiff[is0]) = matDiff;       
1678       }
1679     }
1680     for (Int_t is0=0;is0<72; is0++) {
1681       if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]);
1682       if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]);
1683        //
1684       (*cstream)<<"GAlign"<<
1685         "iter="<<iter<<
1686         "s0="<<is0<<
1687         "m6.="<<arrayAlign[is0]<<
1688         "\n";
1689     }    
1690   }
1691   delete cstream;
1692   for (Int_t isec=0;isec<72;isec++){
1693     fCombinedMatrixArray6.AddAt(arrayAlign[isec],isec);
1694     delete arrayAlignDiff[isec];
1695   }
1696 }
1697
1698
1699 /*
1700   
1701
1702 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1703 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1704 AliXRDPROOFtoolkit tool;
1705 TChain * chainTr = tool.MakeChain("align.txt","Track",0,10200);
1706 chainTr->Lookup();
1707 TChain * chainTracklet = tool.MakeChain("align.txt","Tracklet",0,20);
1708 chainTracklet->Lookup();
1709
1710
1711 TCut cutS0("sqrt(tp2.fC[0]+tp2.fC[0])<0.6");
1712 TCut cutS1("sqrt(tp2.fC[2]+tp2.fC[2])<0.6");
1713 TCut cutS2("sqrt(tp2.fC[5]+tp2.fC[5])<0.04");
1714 TCut cutS3("sqrt(tp2.fC[9]+tp2.fC[9])<0.02");
1715 TCut cutS4("sqrt(tp2.fC[14]+tp2.fC[14])<0.5");
1716 // resolution cuts
1717 TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
1718
1719 TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1720 TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1721 TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.02");
1722 TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.01");
1723 TCut cutE("abs(tp2.fP[1])<235");
1724 TCut cutP=cutP0+cutP1+cutP2+cutP3+cutE;
1725
1726
1727 //
1728 //
1729 TCut cutA =  cutP+cutS;
1730 chainTr->Draw(">>listEL",cutA,"entryList");
1731 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1732 chainTr->SetEntryList(elist);
1733
1734
1735
1736
1737
1738
1739 TCut cutS("s1%36==s2%36");
1740
1741 TCut cutN("c1>32&&c2>60");
1742 TCut cutC0("sqrt(tp2.fC[0])<1");
1743
1744 TCut cutX("abs(tp2.fX-133.6)<2");
1745
1746 TCut cutA = cutP+cutN;
1747
1748
1749 TCut cutY("abs(vcY.fElements-vtY.fElements)<0.3&&vcY.fElements!=0")
1750 TCut cutZ("abs(vcZ.fElements-vtZ.fElements)<0.3&&vcZ.fElements!=0")
1751
1752
1753
1754
1755
1756
1757 TCut cutA =  cutP+cutS;
1758 chainTracklet->Draw(">>listEL",cutA,"entryList");
1759 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1760 chainTracklet->SetEntryList(elist);
1761
1762 //
1763 TVectorD * vec1 = 0;
1764 TVectorD * vec2 = 0;
1765 AliExternalTrackParam * tp1 = 0;
1766 AliExternalTrackParam * tp2 = 0;
1767
1768 Int_t      s1 = 0;
1769 Int_t      s2 = 0;
1770 chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1771 chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1772 chainTracklet->GetBranch("s1")->SetAddress(&s1);
1773 chainTracklet->GetBranch("s2")->SetAddress(&s2);
1774
1775
1776 AliTPCcalibAlign align;
1777 {
1778 for (Int_t i=0; i< elist->GetN(); i++){
1779 //for (Int_t i=0; i<100000; i++){
1780 chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
1781 chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
1782 chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1783 chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1784 chainTracklet->GetBranch("s1")->SetAddress(&s1);
1785 chainTracklet->GetBranch("s2")->SetAddress(&s2);
1786
1787 chainTracklet->GetEntry(i);
1788 if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i,tentry, s1,s2);
1789 //vec1.Print();
1790 TLinearFitter * fitter = align.GetOrMakeFitter6(s1,s2);
1791 if (fitter) align.Process6(vec1->GetMatrixArray(),vec2->GetMatrixArray(),fitter);
1792 }
1793 }
1794
1795 */