During simulation: fill STU region w/ non null time sums
[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 "AliTPCROC.h"
109 #include "AliTPCPointCorrection.h"
110 #include "AliTrackPointArray.h"
111
112 #include "AliExternalTrackParam.h"
113 #include "AliESDEvent.h"
114 #include "AliESDfriend.h"
115 #include "AliESDtrack.h"
116
117 #include "AliTPCTracklet.h"
118 #include "TH1D.h"
119 #include "TH2F.h"
120 #include "THnSparse.h"
121 #include "TVectorD.h"
122 #include "TTreeStream.h"
123 #include "TFile.h"
124 #include "TTree.h"
125 #include "TF1.h"
126 #include "TGraphErrors.h"
127 #include "AliTPCclusterMI.h"
128 #include "AliTPCseed.h"
129 #include "AliTracker.h"
130 #include "TClonesArray.h"
131 #include "AliLog.h"
132 #include "TFile.h"
133 #include "TProfile.h"
134 #include "TCanvas.h"
135 #include "TDatabasePDG.h"
136
137 #include "TTreeStream.h"
138 #include "Riostream.h"
139 #include <sstream>
140 using namespace std;
141
142 AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0;
143 ClassImp(AliTPCcalibAlign)
144
145
146
147
148 AliTPCcalibAlign* AliTPCcalibAlign::Instance()
149 {
150   //
151   // Singleton implementation
152   // Returns an instance of this class, it is created if neccessary
153   //
154   if (fgInstance == 0){
155     fgInstance = new AliTPCcalibAlign();
156   }
157   return fgInstance;
158 }
159
160
161
162
163 AliTPCcalibAlign::AliTPCcalibAlign()
164   :  AliTPCcalibBase(),
165      fDphiHistArray(72*72),
166      fDthetaHistArray(72*72),
167      fDyHistArray(72*72),
168      fDzHistArray(72*72),
169      //
170      fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
171      fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
172      fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
173      fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
174      fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
175      fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     
176      fFitterArray12(72*72),
177      fFitterArray9(72*72),
178      fFitterArray6(72*72),
179      fMatrixArray12(72*72),
180      fMatrixArray9(72*72),
181      fMatrixArray6(72*72),
182      fCombinedMatrixArray6(72),
183      fNoField(kFALSE),
184      fXIO(0),
185      fXmiddle(0),
186      fXquadrant(0),
187      fArraySectorIntParam(36), // array of sector alignment parameters
188      fArraySectorIntCovar(36), // array of sector alignment covariances 
189      //
190      // Kalman filter for global alignment
191      //
192      fSectorParamA(0),     // Kalman parameter   for A side
193      fSectorCovarA(0),     // Kalman covariance  for A side 
194      fSectorParamC(0),     // Kalman parameter   for A side
195      fSectorCovarC(0),     // Kalman covariance  for A side 
196      fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
197 {
198   //
199   // Constructor
200   //
201   for (Int_t i=0;i<72*72;++i) {
202     fPoints[i]=0;
203   }
204   AliTPCROC * roc = AliTPCROC::Instance();
205   fXquadrant = roc->GetPadRowRadii(36,53);
206   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
207   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
208   fClusterDelta[0]=0;   // cluster residuals -  Y
209   fClusterDelta[1]=0;   // cluster residuals -  Z
210   
211   
212   fTrackletDelta[0]=0;   // tracklet residuals
213   fTrackletDelta[1]=0;   // tracklet residuals
214   fTrackletDelta[2]=0;   // tracklet residuals 
215   fTrackletDelta[3]=0;   // tracklet residuals
216 }
217
218 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
219   :AliTPCcalibBase(),  
220    fDphiHistArray(72*72),
221    fDthetaHistArray(72*72),
222    fDyHistArray(72*72),
223    fDzHistArray(72*72),
224    fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
225    fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
226    fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
227    fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
228    fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
229    fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     //
230    fFitterArray12(72*72),
231    fFitterArray9(72*72),
232    fFitterArray6(72*72),
233    fMatrixArray12(72*72),
234    fMatrixArray9(72*72),
235    fMatrixArray6(72*72),
236    fCombinedMatrixArray6(72),
237    fNoField(kFALSE),
238    fXIO(0),
239    fXmiddle(0),
240    fXquadrant(0),
241    fArraySectorIntParam(36), // array of sector alignment parameters
242    fArraySectorIntCovar(36), // array of sector alignment covariances 
243    //
244    // Kalman filter for global alignment
245    //
246    fSectorParamA(0),     // Kalman parameter   for A side
247    fSectorCovarA(0),     // Kalman covariance  for A side 
248    fSectorParamC(0),     // Kalman parameter   for A side
249    fSectorCovarC(0),     // Kalman covariance  for A side      
250    fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
251
252 {
253   //
254   // Constructor
255   //
256   SetName(name);
257   SetTitle(title);
258   for (Int_t i=0;i<72*72;++i) {
259     fPoints[i]=0;
260   }
261   AliTPCROC * roc = AliTPCROC::Instance();
262   fXquadrant = roc->GetPadRowRadii(36,53);
263   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
264   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
265   fClusterDelta[0]=0;   // cluster residuals
266   fClusterDelta[1]=0;   // cluster residuals
267
268   fTrackletDelta[0]=0;   // tracklet residuals
269   fTrackletDelta[1]=0;   // tracklet residuals
270   fTrackletDelta[2]=0;   // tracklet residuals 
271   fTrackletDelta[3]=0;   // tracklet residuals
272 }
273
274
275 AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
276   :AliTPCcalibBase(align),  
277    fDphiHistArray(align.fDphiHistArray),
278    fDthetaHistArray(align.fDthetaHistArray),
279    fDyHistArray(align.fDyHistArray),
280    fDzHistArray(align.fDzHistArray),
281    fDyPhiHistArray(align.fDyPhiHistArray),      // array of residual histograms  y     -kYPhi
282    fDzThetaHistArray(align.fDzThetaHistArray),    // array of residual histograms  z-z   -kZTheta
283    fDphiZHistArray(align.fDphiZHistArray),      // array of residual histograms  phi   -kPhiz
284    fDthetaZHistArray(align.fDthetaZHistArray),    // array of residual histograms  theta -kThetaz
285    fDyZHistArray(align.fDyZHistArray),        // array of residual histograms  y     -kYz
286    fDzZHistArray(align.fDzZHistArray),        // array of residual histograms  z     -kZz     
287    //
288    fFitterArray12(align.fFitterArray12),
289    fFitterArray9(align.fFitterArray9),
290    fFitterArray6(align.fFitterArray6),
291    
292    fMatrixArray12(align.fMatrixArray12),
293    fMatrixArray9(align.fMatrixArray9),
294    fMatrixArray6(align.fMatrixArray6),
295    fCombinedMatrixArray6(align.fCombinedMatrixArray6),
296    fNoField(align.fNoField),
297    fXIO(align.fXIO),   
298    fXmiddle(align.fXmiddle),   
299    fXquadrant(align.fXquadrant),   
300    fArraySectorIntParam(align.fArraySectorIntParam), // array of sector alignment parameters
301    fArraySectorIntCovar(align.fArraySectorIntCovar), // array of sector alignment covariances 
302    fSectorParamA(0),     // Kalman parameter   for A side
303    fSectorCovarA(0),     // Kalman covariance  for A side 
304    fSectorParamC(0),     // Kalman parameter   for A side
305    fSectorCovarC(0),      // Kalman covariance  for A side 
306   fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
307   
308 {
309   //
310   // copy constructor - copy also the content
311   //
312   TH1 * his = 0;
313   TObjArray * arr0=0;
314   const TObjArray *arr1=0;
315   for (Int_t index =0; index<72*72; index++){
316     for (Int_t iarray=0;iarray<10; iarray++){
317       if (iarray==kY){
318         arr0 = &fDyHistArray;
319         arr1 = &align.fDyHistArray;
320       }
321       if (iarray==kZ){
322         arr0 = &fDzHistArray;
323         arr1 = &align.fDzHistArray;
324       }
325       if (iarray==kPhi){
326         arr0 = &fDphiHistArray;
327         arr1 = &align.fDphiHistArray;
328       }
329       if (iarray==kTheta){
330         arr0 = &fDthetaHistArray;
331         arr1 = &align.fDthetaHistArray;
332       }
333       if (iarray==kYz){
334         arr0 = &fDyZHistArray;
335         arr1 = &align.fDyZHistArray;
336       }
337       if (iarray==kZz){
338         arr0 = &fDzZHistArray;
339         arr1 = &align.fDzZHistArray;
340       }
341       if (iarray==kPhiZ){
342         arr0 = &fDphiZHistArray;
343         arr1 = &align.fDphiZHistArray;
344       }
345       if (iarray==kThetaZ){
346         arr0 = &fDthetaZHistArray;
347         arr1 = &align.fDthetaZHistArray;
348       }
349
350       if (iarray==kYPhi){
351         arr0 = &fDyPhiHistArray;
352         arr1 = &align.fDyPhiHistArray;
353       }
354       if (iarray==kZTheta){
355         arr0 = &fDzThetaHistArray;
356         arr1 = &align.fDzThetaHistArray;
357       }
358
359       if (arr1->At(index)) {
360         his = (TH1*)arr1->At(index)->Clone();
361         his->SetDirectory(0);
362         arr0->AddAt(his,index);
363       }    
364     }
365   }
366   //
367   //
368   //
369   if (align.fSectorParamA){
370     fSectorParamA = (TMatrixD*)align.fSectorParamA->Clone();
371     fSectorParamA = (TMatrixD*)align.fSectorCovarA->Clone();
372     fSectorParamC = (TMatrixD*)align.fSectorParamA->Clone();
373     fSectorParamC = (TMatrixD*)align.fSectorCovarA->Clone();
374   }
375   fClusterDelta[0]=0;   // cluster residuals
376   fClusterDelta[1]=0;   // cluster residuals
377
378   fTrackletDelta[0]=0;   // tracklet residuals
379   fTrackletDelta[1]=0;   // tracklet residuals
380   fTrackletDelta[2]=0;   // tracklet residuals 
381   fTrackletDelta[3]=0;   // tracklet residuals
382 }
383
384
385 AliTPCcalibAlign::~AliTPCcalibAlign() {
386   //
387   // destructor
388   //
389   fDphiHistArray.SetOwner(kTRUE);    // array of residual histograms  phi      -kPhi
390   fDthetaHistArray.SetOwner(kTRUE);  // array of residual histograms  theta    -kTheta
391   fDyHistArray.SetOwner(kTRUE);      // array of residual histograms  y        -kY
392   fDzHistArray.SetOwner(kTRUE);      // array of residual histograms  z        -kZ
393   //
394   fDyPhiHistArray.SetOwner(kTRUE);      // array of residual histograms  y     -kYPhi
395   fDzThetaHistArray.SetOwner(kTRUE);    // array of residual histograms  z-z   -kZTheta
396   //
397   fDphiZHistArray.SetOwner(kTRUE);      // array of residual histograms  phi   -kPhiz
398   fDthetaZHistArray.SetOwner(kTRUE);    // array of residual histograms  theta -kThetaz
399   fDyZHistArray.SetOwner(kTRUE);        // array of residual histograms  y     -kYz
400   fDzZHistArray.SetOwner(kTRUE);        // array of residual histograms  z     -kZz
401
402   fDphiHistArray.Delete();    // array of residual histograms  phi      -kPhi
403   fDthetaHistArray.Delete();  // array of residual histograms  theta    -kTheta
404   fDyHistArray.Delete();      // array of residual histograms  y        -kY
405   fDzHistArray.Delete();      // array of residual histograms  z        -kZ
406   //
407   fDyPhiHistArray.Delete();      // array of residual histograms  y     -kYPhi
408   fDzThetaHistArray.Delete();    // array of residual histograms  z-z   -kZTheta
409   //
410   fDphiZHistArray.Delete();      // array of residual histograms  phi   -kPhiz
411   fDthetaZHistArray.Delete();    // array of residual histograms  theta -kThetaz
412   fDyZHistArray.Delete();        // array of residual histograms  y     -kYz
413   fDzZHistArray.Delete();        // array of residual histograms  z     -kZz
414
415   fFitterArray12.SetOwner(kTRUE);    // array of fitters
416   fFitterArray9.SetOwner(kTRUE);     // array of fitters
417   fFitterArray6.SetOwner(kTRUE);     // array of fitters
418   //
419   fMatrixArray12.SetOwner(kTRUE);    // array of transnformtation matrix
420   fMatrixArray9.SetOwner(kTRUE);     // array of transnformtation matrix
421   fMatrixArray6.SetOwner(kTRUE);     // array of transnformtation matrix 
422   //
423   fFitterArray12.Delete();    // array of fitters
424   fFitterArray9.Delete();     // array of fitters
425   fFitterArray6.Delete();     // array of fitters
426   //
427   fMatrixArray12.Delete();    // array of transnformtation matrix
428   fMatrixArray9.Delete();     // array of transnformtation matrix
429   fMatrixArray6.Delete();     // array of transnformtation matrix 
430
431
432   fArraySectorIntParam.SetOwner(kTRUE); // array of sector alignment parameters
433   fArraySectorIntCovar.SetOwner(kTRUE); // array of sector alignment covariances 
434   fArraySectorIntParam.Delete(); // array of sector alignment parameters
435   fArraySectorIntCovar.Delete(); // array of sector alignment covariances 
436   for (Int_t i=0; i<2; i++){
437     delete fClusterDelta[i];   // cluster residuals
438   }
439
440   for (Int_t i=0; i<4; i++){
441     delete fTrackletDelta[i];   // tracklet residuals
442   }
443
444
445 }
446
447 void AliTPCcalibAlign::Process(AliESDEvent *event) {
448   //
449   // Process pairs of cosmic tracks
450   //
451   if (!fClusterDelta[0])  MakeResidualHistos();
452   if (!fTrackletDelta[0])  MakeResidualHistosTracklet();
453   //
454   fCurrentEvent=event;
455   ExportTrackPoints(event);  // export track points for external calibration 
456   const Int_t kMaxTracks =6;
457   const Int_t kminCl = 40;
458   AliESDfriend *eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
459   if (!eESDfriend) return;
460   Int_t ntracks=event->GetNumberOfTracks(); 
461   Float_t dca0[2];
462   Float_t dca1[2];
463   //
464   //
465   // process seeds
466   //
467   for (Int_t i0=0;i0<ntracks;++i0) {
468     AliESDtrack *track0 = event->GetTrack(i0);
469     AliESDfriendTrack *friendTrack = 0;
470     TObject *calibObject=0;
471     AliTPCseed *seed0 = 0;
472     //
473     friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
474     if (!friendTrack) continue;
475     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
476       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
477     }
478     if (!seed0) continue;
479     fCurrentTrack=track0;
480     fCurrentFriendTrack=friendTrack;
481     fCurrentSeed=seed0;
482     fCurrentEvent=event;
483     ProcessSeed(seed0);
484   }
485   //
486   // process cosmic pairs
487   //
488   if (ntracks>kMaxTracks) return;  
489   //
490   //select pairs - for alignment  
491   for (Int_t i0=0;i0<ntracks;++i0) {
492     AliESDtrack *track0 = event->GetTrack(i0);
493     //    if (track0->GetTPCNcls()<kminCl) continue;
494     track0->GetImpactParameters(dca0[0],dca0[1]);
495     //    if (TMath::Abs(dca0[0])>30) continue;
496     //
497     for (Int_t i1=0;i1<ntracks;++i1) {
498       if (i0==i1) continue;
499       AliESDtrack *track1 = event->GetTrack(i1);
500       //      if (track1->GetTPCNcls()<kminCl) continue;
501       track1->GetImpactParameters(dca1[0],dca1[1]);
502       // fast cuts on dca and theta
503       //      if (TMath::Abs(dca1[0]+dca0[0])>15) continue;
504       //      if (TMath::Abs(dca1[1]-dca0[1])>15) continue;
505       if (TMath::Abs(track0->GetParameter()[3]+track1->GetParameter()[3])>0.1) continue;
506       //
507       AliESDfriendTrack *friendTrack = 0;
508       TObject *calibObject=0;
509       AliTPCseed *seed0 = 0,*seed1=0;
510       //
511       friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
512       if (!friendTrack) continue;
513       for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
514         if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
515       }
516       friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i1);;
517       if (!friendTrack) continue;
518       for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
519         if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
520       }
521       if (!seed0) continue;
522       //
523       //
524       if (!seed1) continue;
525       Int_t nclsectors0[72], nclsectors1[72];
526       for (Int_t isec=0;isec<72;isec++){
527         nclsectors0[isec]=0;
528         nclsectors1[isec]=0;
529       }
530       for (Int_t i=0;i<160;i++){
531         AliTPCclusterMI *c0=seed0->GetClusterPointer(i);
532         AliTPCclusterMI *c1=seed1->GetClusterPointer(i);
533         if (c0)  nclsectors0[c0->GetDetector()]+=1;
534         if (c1)  nclsectors1[c1->GetDetector()]+=1;
535       }
536
537       for (Int_t isec0=0; isec0<72;isec0++){
538         if (nclsectors0[isec0]<kminCl) continue;
539         for (Int_t isec1=0; isec1<72;isec1++){
540           if (nclsectors1[isec1]<kminCl) continue;
541           Int_t s0 = isec0;
542           Int_t s1 = isec1;
543           Double_t parLine0[10];
544           Double_t parLine1[10];
545           TMatrixD par0(4,1),cov0(4,4),par1(4,1),cov1(4,4);
546           Bool_t useInnerOuter = kFALSE;
547           if (s1%36!=s0%36) useInnerOuter = fUseInnerOuter;  // for left - right alignment both sectors refit can be used if specified
548           Int_t nl0 = RefitLinear(seed0,s0, parLine0, s0,par0,cov0,fXIO,useInnerOuter);
549           Int_t nl1 = RefitLinear(seed1,s1, parLine1, s0,par1,cov1,fXIO,useInnerOuter);
550           parLine0[0]=0;  // reference frame in IO boundary
551           parLine1[0]=0;
552           //      if (nl0<kminCl || nl1<kminCl) continue;
553           //
554           //
555           Bool_t isOK=kTRUE;
556           if (TMath::Min(nl0,nl1)<kminCl) isOK=kFALSE;
557           // apply selection criteria
558           //
559           Float_t dp0,dp1,dp3;
560           Float_t pp0,pp1,pp3;
561           dp0=par0(0,0)-par1(0,0); 
562           dp1=par0(1,0)-par1(1,0); 
563           dp3=par0(3,0)-par1(3,0); 
564           pp0=dp0/TMath::Sqrt(cov0(0,0)+cov1(0,0)+0.1*0.1);
565           pp1=dp1/TMath::Sqrt(cov0(1,1)+cov1(1,1)+0.0015*0.0015);
566           pp3=dp3/TMath::Sqrt(cov0(3,3)+cov1(3,3)+0.0015*0.0015);
567           //
568           if (TMath::Abs(dp0)>1.0)  isOK=kFALSE;
569           if (TMath::Abs(dp1)>0.02) isOK=kFALSE;
570           if (TMath::Abs(dp3)>0.02) isOK=kFALSE;
571           if (TMath::Abs(pp0)>6)  isOK=kFALSE;
572           if (TMath::Abs(pp1)>6) isOK=kFALSE;
573           if (TMath::Abs(pp3)>6) isOK=kFALSE;     
574           //
575           if (isOK){
576             FillHisto(parLine0,parLine1,s0,s1);  
577             ProcessAlign(parLine0,parLine1,s0,s1);
578             UpdateKalman(s0,s1,par0, cov0, par1, cov1);
579           }
580           if (fStreamLevel>0){
581             TTreeSRedirector *cstream = GetDebugStreamer();
582             if (cstream){
583               (*cstream)<<"cosmic"<<
584                 "isOK="<<isOK<<
585                 "s0="<<s0<<
586                 "s1="<<s1<<
587                 "nl0="<<nl0<<
588                 "nl1="<<nl1<<
589                 "p0.="<<&par0<<
590                 "p1.="<<&par1<<
591                 "c0.="<<&cov0<<
592                 "c1.="<<&cov1<<
593                 "\n";
594             }
595           }
596         }
597       }
598     }
599   }
600 }
601
602 void  AliTPCcalibAlign::ExportTrackPoints(AliESDEvent *event){
603   //
604   // Export track points for alignment - calibration
605   // export space points for pairs of tracks if possible
606   //
607   AliESDfriend *eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
608   if (!eESDfriend) return;
609   Int_t ntracks=event->GetNumberOfTracks();
610   Int_t kMaxTracks=4;   // maximal number of tracks for cosmic pairs
611   Int_t kMinVertexTracks=5;   // maximal number of tracks for vertex mesurement
612
613   //cuts
614   const Int_t kminCl     = 60;
615   const Int_t kminClSum  = 120;
616   //const Double_t kDistY  = 5;
617   // const Double_t kDistZ  = 40;
618   const Double_t kDistTh = 0.05;
619   const Double_t kDistThS = 0.002;
620   const Double_t kDist1Pt = 0.1;
621   const Double_t kMaxD0 =3;  // max distance to the primary vertex
622   const Double_t kMaxD1 =5;  // max distance to the primary vertex
623   const AliESDVertex *tpcVertex = 0;
624   // get the primary vertex TPC
625   if (ntracks>kMinVertexTracks) {
626     tpcVertex = event->GetPrimaryVertexSPD();
627     if (tpcVertex->GetNContributors()<kMinVertexTracks) tpcVertex=0;
628   }
629   //
630   Float_t dca0[2];
631   //  Float_t dca1[2];
632   Int_t index0=0,index1=0;
633   //
634   for (Int_t i0=0;i0<ntracks;++i0) {
635     AliESDtrack *track0 = event->GetTrack(i0);
636     if (!track0) continue;
637     if ((track0->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
638     if (track0->GetOuterParam()==0) continue;
639     if (track0->GetInnerParam()==0) continue;
640     if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt()-track0->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
641     if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
642     if (TMath::Abs(track0->GetInnerParam()->GetTgl()-track0->GetOuterParam()->GetTgl())>kDistThS) continue;
643     AliESDtrack *track1P = 0;
644     if (track0->GetTPCNcls()<kminCl) continue;
645     track0->GetImpactParameters(dca0[0],dca0[1]);
646     index0=i0;
647     index1=-1;
648     //
649     if (ntracks<kMaxTracks) for (Int_t i1=i0+1;i1<ntracks;++i1) {
650       if (i0==i1) continue;
651       AliESDtrack *track1 = event->GetTrack(i1);
652       if (!track1) continue;
653       if ((track1->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
654       if (track1->GetOuterParam()==0) continue;
655       if (track1->GetInnerParam()==0) continue;
656       if (track1->GetTPCNcls()<kminCl) continue;
657       if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt()-track1->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
658       if (TMath::Abs(track1->GetInnerParam()->GetTgl()-track1->GetOuterParam()->GetTgl())>kDistThS) continue;
659       if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
660       //track1->GetImpactParameters(dca1[0],dca1[1]);
661       //if (TMath::Abs(dca1[0]-dca0[0])>kDistY) continue;
662       //if (TMath::Abs(dca1[1]-dca0[1])>kDistZ) continue;
663       if (TMath::Abs(track0->GetTgl()+track1->GetTgl())>kDistTh) continue;
664       if (TMath::Abs(track0->GetSigned1Pt()+track1->GetSigned1Pt())>kDist1Pt) continue;
665       track1P = track1;
666       index1=i1;
667     }
668     AliESDfriendTrack *friendTrack = 0;
669     TObject *calibObject=0;
670     AliTPCseed *seed0 = 0,*seed1=0;
671     //
672     friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index0);;
673     if (!friendTrack) continue;
674     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
675       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
676     }
677     if (index1>0){
678       friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index1);;
679       if (!friendTrack) continue;
680       for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
681         if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
682       }
683     }
684     //
685     Int_t npoints=0, ncont=0;
686     if (seed0) {npoints+=seed0->GetNumberOfClusters(); ncont++;}
687     if (seed1) {npoints+=seed1->GetNumberOfClusters(); ncont++;}
688     if (npoints<kminClSum) continue;    
689     Int_t cpoint=0;
690     AliTrackPointArray array(npoints);    
691     if (tpcVertex){
692       Double_t dxyz[3]={0,0,0};
693       Double_t dc[6]={0,0,0};
694       tpcVertex->GetXYZ(dxyz);
695       tpcVertex->GetCovarianceMatrix(dc);
696       Float_t xyz[3]={dxyz[0],dxyz[1],dxyz[2]};
697       Float_t cov[6]={dc[0]+1,dc[1],dc[2]+1,dc[3],dc[4],dc[5]+100.}; 
698       AliTrackPoint point(xyz,cov,73);  // add point to not existing volume
699       Float_t dtpc[2],dcov[3];
700       track0->GetImpactParametersTPC(dtpc,dcov);
701       if (TMath::Abs(dtpc[0])>kMaxD0) continue;
702       if (TMath::Abs(dtpc[1])>kMaxD1) continue;
703     }
704
705     if (seed0) for (Int_t icl = 0; icl<160; icl++){
706       AliTPCclusterMI *cluster=seed0->GetClusterPointer(icl);
707       if (!cluster) continue;
708       Float_t xyz[3];
709       Float_t cov[6];
710       cluster->GetGlobalXYZ(xyz);
711       cluster->GetGlobalCov(cov);
712       AliTrackPoint point(xyz,cov,cluster->GetDetector());
713       array.AddPoint(npoints, &point); 
714       if (cpoint>=npoints) continue;  //shoul not happen
715       array.AddPoint(cpoint, &point);
716       cpoint++;
717     }
718     if (seed1) for (Int_t icl = 0; icl<160; icl++){
719       AliTPCclusterMI *cluster=seed1->GetClusterPointer(icl);
720       if (!cluster) continue;
721       Float_t xyz[3];
722       Float_t cov[6];
723       cluster->GetGlobalXYZ(xyz);
724       cluster->GetGlobalCov(cov);
725       AliTrackPoint point(xyz,cov,cluster->GetDetector());
726       array.AddPoint(npoints, &point);
727       if (cpoint>=npoints) continue;  //shoul not happen
728       array.AddPoint(cpoint, &point);
729       cpoint++;
730     }
731     //
732     //
733     //
734     TTreeSRedirector *cstream = GetDebugStreamer();
735     if (cstream){
736       Bool_t isVertex=(tpcVertex)? kTRUE:kFALSE;
737       Double_t tof0=track0->GetTOFsignal();
738       Double_t tof1=(track1P) ?  track1P->GetTOFsignal(): 0;
739       static AliExternalTrackParam dummy;
740       AliExternalTrackParam *p0In  = &dummy;
741       AliExternalTrackParam *p1In  = &dummy;
742       AliExternalTrackParam *p0Out = &dummy;
743       AliExternalTrackParam *p1Out = &dummy;
744       AliESDVertex vdummy;
745       AliESDVertex *pvertex= (tpcVertex)? (AliESDVertex *)tpcVertex: &vdummy;
746       if (track0) {
747         p0In= new AliExternalTrackParam(*track0);
748         p0Out=new AliExternalTrackParam(*(track0->GetOuterParam()));
749       }
750       if (track1P) {
751         p1In= new AliExternalTrackParam(*track1P);
752         p1Out=new AliExternalTrackParam(*(track1P->GetOuterParam()));
753       }
754
755       (*cstream)<<"trackPoints"<<
756         "run="<<fRun<<              //  run number
757         "event="<<fEvent<<          //  event number
758         "time="<<fTime<<            //  time stamp of event
759         "trigger="<<fTrigger<<      //  trigger
760         "triggerClass="<<&fTriggerClass<<      //  trigger
761         "mag="<<fMagF<<             //  magnetic field
762         "pvertex.="<<pvertex<<      // vertex
763         //
764         "isVertex="<<isVertex<<     // flag is with prim vertex
765         "tof0="<<tof0<<             // tof signal 0
766         "tof1="<<tof1<<             // tof signal 1
767         "seed0.="<<seed0<<          // track info
768         "ntracks="<<ntracks<<       // number of tracks
769         "ncont="<<ncont<<           // number of contributors
770         "p0In.="<<p0In<<              // track parameters0 
771         "p1In.="<<p1In<<              // track parameters1
772         "p0Out.="<<p0Out<<          // track parameters0
773         "p1Out.="<<p1Out<<          // track parameters0
774         "p.="<<&array<<
775         "\n";
776     }
777   }  
778 }
779
780
781
782
783 void AliTPCcalibAlign::ProcessSeed(AliTPCseed *seed) {
784   //
785   // 
786   //
787   // make a kalman tracklets out of seed
788   //
789   UpdateClusterDeltaField(seed);
790   TObjArray tracklets=
791     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
792                                     kFALSE,20,4);
793   tracklets.SetOwner();
794   Int_t ntracklets = tracklets.GetEntries();
795   if (ntracklets<2) return;
796   //
797   //
798   for (Int_t i1=0;i1<ntracklets;i1++)
799     for (Int_t i2=0;i2<ntracklets;i2++){
800       if (i1==i2) continue;
801       AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[i1]);
802       AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[i2]);
803       AliExternalTrackParam *common1=0,*common2=0;
804       if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2)){
805         ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());      
806         UpdateAlignSector(seed,t1->GetSector());
807       }
808       delete common1;
809       delete common2;
810     }  
811 }
812
813 void AliTPCcalibAlign::Analyze(){
814   //
815   // Analyze function 
816   //
817   EvalFitters();
818 }
819
820
821 void AliTPCcalibAlign::Terminate(){
822   //
823   // Terminate function
824   // call base terminate + Eval of fitters
825   //
826   Info("AliTPCcalibAlign","Terminate");
827   EvalFitters();
828   AliTPCcalibBase::Terminate();
829 }
830
831
832 void AliTPCcalibAlign::UpdatePointCorrection(AliTPCPointCorrection * correction){
833   //
834   // Update point correction with alignment coefficients
835   //
836   for (Int_t isec=0;isec<36;isec++){
837     TMatrixD * matCorr = (TMatrixD*)(correction->fArraySectorIntParam.At(isec));
838     TMatrixD * matAlign = (TMatrixD*)(fArraySectorIntParam.At(isec));
839     TMatrixD * matAlignCovar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
840     if (!matAlign) continue;
841     if (!matCorr) {
842       correction->fArraySectorIntParam.AddAt(matAlign->Clone(),isec); 
843       correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec); 
844       continue;
845     }
846     (*matCorr)+=(*matAlign);
847     correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec); 
848   }
849   //
850
851 }
852
853
854 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
855                                         const AliExternalTrackParam &tp2,
856                                         const AliTPCseed * seed,
857                                         Int_t s1,Int_t s2) {
858   //
859   // Process function to fill fitters
860   //
861   if (!seed) return;
862   Double_t t1[10],t2[10];
863   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
864   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
865   x1   =tp1.GetX();
866   y1   =tp1.GetY();
867   z1   =tp1.GetZ();
868   Double_t snp1=tp1.GetSnp();
869   dydx1=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
870   Double_t tgl1=tp1.GetTgl();
871   // dz/dx = 1/(cos(theta)*cos(phi))
872   dzdx1=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
873   x2   =tp2.GetX();
874   y2   =tp2.GetY();
875   z2   =tp2.GetZ();
876   Double_t snp2=tp2.GetSnp();
877   dydx2=snp2/TMath::Sqrt((1.-snp2)*(1.+snp2));
878   Double_t tgl2=tp2.GetTgl();
879   dzdx2=tgl2/TMath::Sqrt((1.-snp2)*(1.+snp2));
880   
881   //
882   // Kalman parameters
883   //
884   t1[0]-=fXIO;
885   t2[0]-=fXIO;
886   // errors
887   t1[5]=0; t2[5]=0;
888   t1[6]=TMath::Sqrt(tp1.GetSigmaY2());
889   t1[7]=TMath::Sqrt(tp1.GetSigmaSnp2());
890   t1[8]=TMath::Sqrt(tp1.GetSigmaZ2()); 
891   t1[9]=TMath::Sqrt(tp1.GetSigmaTgl2());
892   
893   t2[6]=TMath::Sqrt(tp2.GetSigmaY2());
894   t2[7]=TMath::Sqrt(tp2.GetSigmaSnp2()); 
895   t2[8]=TMath::Sqrt(tp2.GetSigmaZ2());
896   t2[9]=TMath::Sqrt(tp2.GetSigmaTgl2());
897   //
898   // linear parameters
899   //
900   Double_t parLine1[10];
901   Double_t parLine2[10];
902   TMatrixD par1(4,1),cov1(4,4),par2(4,1),cov2(4,4);
903   Bool_t useInnerOuter = kFALSE;
904   if (s1%36!=s2%36) useInnerOuter = fUseInnerOuter;  // for left - right alignment bot sectors refit can be used if specified
905   Int_t nl1 = RefitLinear(seed,s1, parLine1, s1,par1,cov1,tp1.GetX(), useInnerOuter);
906   Int_t nl2 = RefitLinear(seed,s2, parLine2, s1,par2,cov2,tp1.GetX(), useInnerOuter);
907   parLine1[0]=tp1.GetX()-fXIO;   // parameters in  IROC-OROC boundary
908   parLine2[0]=tp1.GetX()-fXIO;   // parameters in  IROC-OROC boundary
909   //
910   //
911   //
912   Int_t accept       =   AcceptTracklet(tp1,tp2);  
913   Int_t acceptLinear =   AcceptTracklet(parLine1,parLine2);
914
915
916   if (fStreamLevel>1){
917     TTreeSRedirector *cstream = GetDebugStreamer();
918     if (cstream){
919       static TVectorD vec1(5);
920       static TVectorD vec2(5);
921       static TVectorD vecL1(9);
922       static TVectorD vecL2(9);
923       vec1.SetElements(t1);
924       vec2.SetElements(t2);
925       vecL1.SetElements(parLine1);
926       vecL2.SetElements(parLine2);
927       AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
928       AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
929       (*cstream)<<"Tracklet"<<
930         "accept="<<accept<<
931         "acceptLinear="<<acceptLinear<<  // accept linear tracklets
932         "run="<<fRun<<              //  run number
933         "event="<<fEvent<<          //  event number
934         "time="<<fTime<<            //  time stamp of event
935         "trigger="<<fTrigger<<      //  trigger
936         "triggerClass="<<&fTriggerClass<<      //  trigger
937         "mag="<<fMagF<<             //  magnetic field
938         "isOK="<<accept<<           //  flag - used for alignment
939         "tp1.="<<p1<<
940         "tp2.="<<p2<<
941         "v1.="<<&vec1<<
942         "v2.="<<&vec2<<
943         "s1="<<s1<<
944         "s2="<<s2<<
945         "nl1="<<nl1<<       // linear fit - n points
946         "nl2="<<nl2<<       // linear fit - n points
947         "vl1.="<<&vecL1<<   // linear fits
948         "vl2.="<<&vecL2<<   // linear fits
949         "\n";
950     }
951   }
952   if (TMath::Abs(fMagF)<0.005){
953     //
954     // use Linear fit
955     //
956     if (nl1>10 && nl2>10 &&(acceptLinear==0)){
957       ProcessDiff(tp1,tp2, seed,s1,s2);
958       if (TMath::Abs(parLine1[2])<0.8 &&TMath::Abs(parLine1[2])<0.8 ){ //angular cut
959         FillHisto(parLine1,parLine2,s1,s2);  
960         ProcessAlign(parLine1,parLine2,s1,s2);
961         FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
962         FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
963         //UpdateKalman(s1,s2,par1, cov1, par2, cov2); - OBSOLETE to be removed - 50 % of time here
964       }
965     }
966   }
967   if (accept>0) return;
968   //
969   // fill resolution histograms - previous cut included
970   if (TMath::Abs(fMagF)>0.005){
971     //
972     // use Kalman if mag field
973     //
974     ProcessDiff(tp1,tp2, seed,s1,s2);
975     FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
976     FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
977     FillHisto(t1,t2,s1,s2);  
978     ProcessAlign(t1,t2,s1,s2);
979   }
980 }
981
982 void AliTPCcalibAlign::ProcessAlign(Double_t * t1,
983                                     Double_t * t2,
984                                     Int_t s1,Int_t s2){
985   //
986   // Do intersector alignment
987   //
988   //Process12(t1,t2,GetOrMakeFitter12(s1,s2));
989   //Process9(t1,t2,GetOrMakeFitter9(s1,s2));
990   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
991   ++fPoints[GetIndex(s1,s2)];
992 }
993
994
995
996 Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
997                                        const AliExternalTrackParam &p2) const
998 {
999
1000   //
1001   // Accept pair of tracklets?
1002   //
1003   /*
1004   // resolution cuts
1005   TCut cutS0("sqrt(tp2.fC[0]+tp1.fC[0])<0.2");
1006   TCut cutS1("sqrt(tp2.fC[2]+tp1.fC[2])<0.2");
1007   TCut cutS2("sqrt(tp2.fC[5]+tp1.fC[5])<0.01");
1008   TCut cutS3("sqrt(tp2.fC[9]+tp1.fC[9])<0.01");
1009   TCut cutS4("sqrt(tp2.fC[14]+tp1.fC[14])<0.25");
1010   TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
1011   //
1012   // parameters matching cuts
1013   TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1014   TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1015   TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03");
1016   TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03");
1017   TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
1018   TCut cutPP4("abs(tp1.fP[4]-tp2.fP[4])/sqrt(tp2.fC[14]+tp1.fC[14])<3");
1019   TCut cutP=cutP0+cutP1+cutP2+cutP3+cutP4+cutPP4;
1020   */  
1021   //
1022   // resolution cuts
1023   Int_t reject=0;
1024   const Double_t *cp1 = p1.GetCovariance();
1025   const Double_t *cp2 = p2.GetCovariance();
1026   if (TMath::Sqrt(cp1[0]+cp2[0])>0.2)  reject|=1;;
1027   if (TMath::Sqrt(cp1[2]+cp2[2])>0.2)  reject|=2;
1028   if (TMath::Sqrt(cp1[5]+cp2[5])>0.01) reject|=4;
1029   if (TMath::Sqrt(cp1[9]+cp2[9])>0.01) reject|=8;
1030   if (TMath::Sqrt(cp1[14]+cp2[14])>0.2) reject|=16;
1031
1032   //parameters difference
1033   const Double_t *tp1 = p1.GetParameter();
1034   const Double_t *tp2 = p2.GetParameter();
1035   if (TMath::Abs(tp1[0]-tp2[0])>0.6) reject|=32;
1036   if (TMath::Abs(tp1[1]-tp2[1])>0.6) reject|=64;
1037   if (TMath::Abs(tp1[2]-tp2[2])>0.03) reject|=128;
1038   if (TMath::Abs(tp1[3]-tp2[3])>0.03) reject|=526;
1039   if (TMath::Abs(tp1[4]-tp2[4])>0.4) reject|=1024;
1040   if (TMath::Abs(tp1[4]-tp2[4])/TMath::Sqrt(cp1[14]+cp2[14])>4) reject|=2048;
1041   
1042   //
1043   if (TMath::Abs(tp2[1])>235) reject|=2*4096;
1044   
1045   if (fNoField){
1046     
1047   }
1048
1049   return reject;
1050 }
1051
1052
1053 Int_t  AliTPCcalibAlign::AcceptTracklet(const Double_t *t1, const Double_t *t2) const
1054 {
1055   //
1056   // accept tracklet  - 
1057   //  dist cut + 6 sigma cut 
1058   //
1059   Double_t dy     = t2[1]-t1[1];
1060   Double_t dphi   = t2[2]-t1[2];
1061   Double_t dz     = t2[3]-t1[3];
1062   Double_t dtheta = t2[4]-t1[4];
1063   //
1064   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]+0.05*0.05);
1065   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]+0.001*0.001);
1066   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]+0.05*0.05);
1067   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]+0.001*0.001);
1068   //
1069   Int_t reject=0;
1070   if (TMath::Abs(dy)>1.)         reject|=2;
1071   if (TMath::Abs(dphi)>0.1)      reject|=4;
1072   if (TMath::Abs(dz)>1.)         reject|=8;
1073   if (TMath::Abs(dtheta)>0.1)    reject|=16;
1074   //
1075   if (TMath::Abs(dy/sy)>6)         reject|=32;
1076   if (TMath::Abs(dphi/sdydx)>6)    reject|=64;
1077   if (TMath::Abs(dz/sz)>6)         reject|=128;
1078   if (TMath::Abs(dtheta/sdzdx)>6)  reject|=256;
1079   return reject;
1080 }
1081
1082
1083 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
1084                                     const AliExternalTrackParam &t2,
1085                                     const AliTPCseed *seed,
1086                                     Int_t s1,Int_t s2)
1087 {
1088   //
1089   // Process local residuals function
1090   // 
1091   TVectorD vecX(160);
1092   TVectorD vecY(160);
1093   TVectorD vecZ(160);
1094   TVectorD vecClY(160);
1095   TVectorD vecClZ(160);
1096   TClonesArray arrCl("AliTPCclusterMI",160);
1097   arrCl.ExpandCreateFast(160);
1098   Int_t count1=0, count2=0;
1099   
1100   for (Int_t i=0;i<160;++i) {
1101     AliTPCclusterMI *c=seed->GetClusterPointer(i);
1102     vecX[i]=0;
1103     vecY[i]=0;
1104     vecZ[i]=0;
1105     if (!c) continue;
1106     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
1107     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
1108     vecClY[i] = c->GetY();
1109     vecClZ[i] = c->GetZ();
1110     cl=*c;
1111     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
1112     if (c->GetDetector()==s1) ++count1;
1113     if (c->GetDetector()==s2) ++count2;
1114     Double_t gxyz[3],xyz[3];
1115     t1.GetXYZ(gxyz);
1116     Float_t bz = AliTracker::GetBz(gxyz);
1117     par->GetYAt(c->GetX(), bz, xyz[1]);
1118     par->GetZAt(c->GetX(), bz, xyz[2]);
1119     vecX[i] = c->GetX();
1120     vecY[i]= xyz[1];
1121     vecZ[i]= xyz[2];
1122   }
1123   //
1124   //
1125   if (fStreamLevel>5 && count1>10 && count2>10){
1126     //
1127     // huge output - cluster residuals to be investigated
1128     //
1129     TTreeSRedirector *cstream = GetDebugStreamer();
1130     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
1131     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
1132     /*
1133       
1134       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");
1135
1136     */
1137
1138     if (cstream){
1139       (*cstream)<<"Track"<<
1140         "run="<<fRun<<              //  run number
1141         "event="<<fEvent<<          //  event number
1142         "time="<<fTime<<            //  time stamp of event
1143         "trigger="<<fTrigger<<      //  trigger
1144         "triggerClass="<<&fTriggerClass<<      //  trigger
1145         "mag="<<fMagF<<             //  magnetic field
1146         "Cl.="<<&arrCl<<
1147         //"tp0.="<<p0<<
1148         "tp1.="<<p1<<
1149         "tp2.="<<p2<<
1150         "vtX.="<<&vecX<<
1151         "vtY.="<<&vecY<<
1152         "vtZ.="<<&vecZ<<
1153         "vcY.="<<&vecClY<<
1154         "vcZ.="<<&vecClZ<<
1155         "s1="<<s1<<
1156         "s2="<<s2<<
1157         "c1="<<count1<<
1158         "c2="<<count2<<
1159         "\n";
1160     }
1161   }
1162 }
1163
1164
1165
1166
1167 void AliTPCcalibAlign::Process12(const Double_t *t1,
1168                                  const Double_t *t2,
1169                                  TLinearFitter *fitter) const
1170 {
1171   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1172   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1173   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1174   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1175   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1176   //
1177   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
1178   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
1179   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
1180
1181
1182
1183   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1184   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1185
1186   //
1187   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1188   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1189   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1190   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1191
1192   Double_t p[12];
1193   Double_t value;
1194
1195   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1196   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
1197   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1198   for (Int_t i=0; i<12;i++) p[i]=0.;
1199   p[3+0] = x1;          // a10
1200   p[3+1] = y1;          // a11
1201   p[3+2] = z1;          // a12
1202   p[9+1] = 1.;          // a13
1203   p[0+1] = y1*dydx2;    // a01
1204   p[0+2] = z1*dydx2;    // a02
1205   p[9+0] = dydx2;       // a03
1206   value  = y2;
1207   fitter->AddPoint(p,value,sy);
1208
1209   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1210   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
1211   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1212   for (Int_t i=0; i<12;i++) p[i]=0.;
1213   p[6+0] = x1;           // a20 
1214   p[6+1] = y1;           // a21
1215   p[6+2] = z1;           // a22
1216   p[9+2] = 1.;           // a23
1217   p[0+1] = y1*dzdx2;     // a01
1218   p[0+2] = z1*dzdx2;     // a02
1219   p[9+0] = dzdx2;        // a03
1220   value  = z2;
1221   fitter->AddPoint(p,value,sz);
1222
1223   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1224   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1225   for (Int_t i=0; i<12;i++) p[i]=0.;
1226   p[3+0] = 1.;           // a10
1227   p[3+1] = dydx1;        // a11
1228   p[3+2] = dzdx1;        // a12
1229   p[0+0] = -dydx2;       // a00
1230   p[0+1] = -dydx1*dydx2; // a01
1231   p[0+2] = -dzdx1*dydx2; // a02
1232   value  = 0.;
1233   fitter->AddPoint(p,value,sdydx);
1234
1235   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1236   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1237   for (Int_t i=0; i<12;i++) p[i]=0.;
1238   p[6+0] = 1;            // a20
1239   p[6+1] = dydx1;        // a21
1240   p[6+2] = dzdx1;        // a22
1241   p[0+0] = -dzdx2;       // a00
1242   p[0+1] = -dydx1*dzdx2; // a01
1243   p[0+2] = -dzdx1*dzdx2; // a02
1244   value  = 0.;
1245   fitter->AddPoint(p,value,sdzdx);
1246 }
1247
1248 void AliTPCcalibAlign::Process9(const Double_t * const t1,
1249                                 const Double_t * const t2,
1250                                 TLinearFitter *fitter) const
1251 {
1252   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1253   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1254   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1255   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1256   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1257   //
1258   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
1259   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
1260   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8] 
1261
1262
1263   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1264   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1265   //
1266   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1267   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1268   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1269   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1270
1271   //
1272   Double_t p[12];
1273   Double_t value;
1274
1275   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1276   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
1277   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1278   for (Int_t i=0; i<12;i++) p[i]=0.;
1279   p[2]   += x1;           // a10
1280   //p[]  +=1;             // a11
1281   p[3]   += z1;           // a12    
1282   p[7]   += 1;            // a13
1283   p[0]   += y1*dydx2;     // a01
1284   p[1]   += z1*dydx2;     // a02
1285   p[6]   += dydx2;        // a03
1286   value   = y2-y1;        //-a11
1287   fitter->AddPoint(p,value,sy);
1288   //
1289   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1290   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
1291   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1292   for (Int_t i=0; i<12;i++) p[i]=0.;
1293   p[4]   += x1;           // a20 
1294   p[5]   += y1;           // a21
1295   //p[]  += z1;           // a22
1296   p[8]   += 1.;           // a23
1297   p[0]   += y1*dzdx2;     // a01
1298   p[1]   += z1*dzdx2;     // a02
1299   p[6]   += dzdx2;        // a03
1300   value  = z2-z1;         //-a22
1301   fitter->AddPoint(p,value,sz);
1302
1303   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1304   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1305   for (Int_t i=0; i<12;i++) p[i]=0.;
1306   p[2]   += 1.;           // a10
1307   //p[]  += dydx1;      // a11
1308   p[3]   += dzdx1;        // a12
1309   //p[]  += -dydx2;       // a00
1310   p[0]   += -dydx1*dydx2; // a01
1311   p[1]   += -dzdx1*dydx2; // a02
1312   value  = -dydx1+dydx2;  // -a11 + a00
1313   fitter->AddPoint(p,value,sdydx);
1314
1315   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1316   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1317   for (Int_t i=0; i<12;i++) p[i]=0.;
1318   p[4]   += 1;            // a20
1319   p[5]   += dydx1;        // a21
1320   //p[]  += dzdx1;        // a22
1321   //p[]  += -dzdx2;       // a00
1322   p[0]   += -dydx1*dzdx2; // a01
1323   p[1]   += -dzdx1*dzdx2; // a02
1324   value  = -dzdx1+dzdx2;  // -a22 + a00
1325   fitter->AddPoint(p,value,sdzdx);
1326 }
1327
1328 void AliTPCcalibAlign::Process6(const Double_t *const t1,
1329                                 const Double_t *const t2,
1330                                 TLinearFitter *fitter) const
1331 {
1332   // x2    =  1  *x1 +-a01*y1 + 0      +a03
1333   // y2    =  a01*x1 + 1  *y1 + 0      +a13
1334   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
1335   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1336   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1337   //
1338   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
1339   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
1340   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5] 
1341
1342   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1343   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1344
1345   //
1346   Double_t sy       = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1347   Double_t sdydx    = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1348   Double_t sz       = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1349   Double_t sdzdx    = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1350
1351  
1352   Double_t p[12];
1353   Double_t value;
1354   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1355   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
1356   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1357   for (Int_t i=0; i<12;i++) p[i]=0.;
1358   p[0]   += x1;           // a10
1359   //p[]  +=1;             // a11
1360   //p[]  += z1;           // a12    
1361   p[4]   += 1;            // a13
1362   p[0]   += -y1*dydx2;    // a01
1363   //p[]  += z1*dydx2;     // a02
1364   p[3]   += dydx2;        // a03
1365   value   = y2-y1;        //-a11
1366   fitter->AddPoint(p,value,sy);
1367   //
1368   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
1369   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
1370   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1371   for (Int_t i=0; i<12;i++) p[i]=0.;
1372   p[1]   += x1;           // a20 
1373   p[2]   += y1;           // a21
1374   //p[]  += z1;           // a22
1375   p[5]   += 1.;           // a23
1376   p[0]   += -y1*dzdx2;    // a01
1377   //p[]   += z1*dzdx2;     // a02
1378   p[3]   += dzdx2;        // a03
1379   value  = z2-z1;         //-a22
1380   fitter->AddPoint(p,value,sz);
1381
1382   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1383   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1384   for (Int_t i=0; i<12;i++) p[i]=0.;
1385   p[0]   += 1.;           // a10
1386   //p[]  += dydx1;      // a11       
1387   //p[]   += dzdx1;        // a12
1388   //p[]  += -dydx2;       // a00
1389   //p[0]   +=  dydx1*dydx2; // a01         FIXME- 0912 MI
1390   //p[]   += -dzdx1*dydx2; // a02
1391   value  = -dydx1+dydx2;  // -a11 + a00
1392   fitter->AddPoint(p,value,sdydx);
1393
1394   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1395   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1396   for (Int_t i=0; i<12;i++) p[i]=0.;
1397   p[1]   += 1;            // a20
1398   //  p[2]   += dydx1;        // a21   FIXME- 0912 MI
1399   //p[]  += dzdx1;        // a22
1400   //p[]  += -dzdx2;       // a00
1401   //p[0]   +=  dydx1*dzdx2; // a01     FIXME- 0912 MI
1402   //p[]  += -dzdx1*dzdx2; // a02
1403   value  = -dzdx1+dzdx2;  // -a22 + a00
1404   fitter->AddPoint(p,value,sdzdx);
1405 }
1406
1407
1408
1409
1410 void AliTPCcalibAlign::EvalFitters(Int_t minPoints) {
1411   //
1412   // Analyze function 
1413   // 
1414   // Perform the fitting using linear fitters
1415   //
1416   TLinearFitter *f;
1417   TFile fff("alignDebug.root","recreate");
1418   for (Int_t s1=0;s1<72;++s1)
1419     for (Int_t s2=0;s2<72;++s2){
1420       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1421         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1422         if (f->Eval()!=0) {
1423           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1424           f->Write(Form("f12_%d_%d",s1,s2));
1425         }else{
1426           f->Write(Form("f12_%d_%d",s1,s2));
1427         }
1428       }
1429       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1430         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1431         if (f->Eval()!=0) {
1432           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1433         }else{
1434           f->Write(Form("f9_%d_%d",s1,s2));
1435         }
1436       }
1437       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1438         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1439         if (f->Eval()!=0) {
1440           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1441         }else{
1442           f->Write(Form("f6_%d_%d",s1,s2));
1443         }
1444       }
1445     }
1446   TMatrixD mat(4,4);
1447   for (Int_t s1=0;s1<72;++s1)
1448     for (Int_t s2=0;s2<72;++s2){
1449       if (GetTransformation12(s1,s2,mat)){
1450         fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2));
1451       }
1452       if (GetTransformation9(s1,s2,mat)){
1453         fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2));
1454       }
1455       if (GetTransformation6(s1,s2,mat)){
1456         fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2));
1457       }
1458     }
1459   //this->Write("align");
1460   
1461 }
1462
1463 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
1464   //
1465   // get or make fitter - general linear transformation
1466   //
1467   static Int_t counter12=0;
1468   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]");
1469   TLinearFitter * fitter = GetFitter12(s1,s2);
1470   if (fitter) return fitter;
1471   //  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]");
1472   fitter =new TLinearFitter(&f12,"");
1473   fitter->StoreData(kFALSE);
1474   fFitterArray12.AddAt(fitter,GetIndex(s1,s2)); 
1475   counter12++;
1476   //  if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
1477   return fitter;
1478 }
1479
1480 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
1481   //
1482   //get or make fitter - general linear transformation - no scaling
1483   // 
1484   static Int_t counter9=0;
1485   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
1486   TLinearFitter * fitter = GetFitter9(s1,s2);
1487   if (fitter) return fitter;
1488   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
1489   fitter =new TLinearFitter(&f9,"");
1490   fitter->StoreData(kFALSE);
1491   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
1492   counter9++;
1493   //  if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
1494   return fitter;
1495 }
1496
1497 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
1498   //
1499   // get or make fitter  - 6 paramater linear tranformation
1500   //                     - no scaling
1501   //                     - rotation x-y
1502   //                     - tilting x-z, y-z
1503   static Int_t counter6=0;
1504   static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1505   TLinearFitter * fitter = GetFitter6(s1,s2);
1506   if (fitter) return fitter;
1507   //  fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1508   fitter=new TLinearFitter(&f6,"");
1509   fitter->StoreData(kFALSE);
1510   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
1511   counter6++;
1512   //  if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
1513   return fitter;
1514 }
1515
1516
1517
1518
1519
1520 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
1521   //
1522   // GetTransformation matrix - 12 paramaters - generael linear transformation
1523   //
1524   if (!GetFitter12(s1,s2))
1525     return false;
1526   else {
1527     TVectorD p(12);
1528     GetFitter12(s1,s2)->GetParameters(p);
1529     a.ResizeTo(4,4);
1530     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
1531     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
1532     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
1533     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1534     return true;
1535   } 
1536 }
1537
1538 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
1539   //
1540   // GetTransformation matrix - 9 paramaters - general linear transformation
1541   //                            No scaling
1542   //
1543   if (!GetFitter9(s1,s2))
1544     return false;
1545   else {
1546     TVectorD p(9);
1547     GetFitter9(s1,s2)->GetParameters(p);
1548     a.ResizeTo(4,4);
1549     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
1550     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
1551     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
1552     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1553     return true;
1554   } 
1555 }
1556
1557 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
1558   //
1559   // GetTransformation matrix - 6  paramaters
1560   //                            3  translation
1561   //                            1  rotation -x-y  
1562   //                            2  tilting x-z y-z
1563   if (!GetFitter6(s1,s2))
1564     return false;
1565   else {
1566     TVectorD p(6);
1567     GetFitter6(s1,s2)->GetParameters(p);
1568     a.ResizeTo(4,4);
1569     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
1570     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
1571     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
1572     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1573     return true;
1574   } 
1575 }
1576
1577 void AliTPCcalibAlign::MakeResidualHistos(){
1578   //
1579   // Make cluster residual histograms
1580   //
1581   Double_t xminTrack[9], xmaxTrack[9];
1582   Int_t    binsTrack[9];
1583   TString  axisName[9],axisTitle[9];
1584   //
1585   // 0 - delta   of interest
1586   // 1 - global  phi in sector number  as float
1587   // 2 - local   x
1588   // 3 - local   ky
1589   // 4 - local   kz
1590   // 
1591   axisName[0]="delta";   axisTitle[0]="#Delta (cm)"; 
1592   if (TMath::Abs(AliTracker::GetBz())<0.01){
1593     binsTrack[0]=60;       xminTrack[0]=-1.2;        xmaxTrack[0]=1.2; 
1594   }else{
1595     binsTrack[0]=60;       xminTrack[0]=-0.6;        xmaxTrack[0]=0.6; 
1596   }
1597   //
1598   axisName[1]="sector";   axisTitle[1]="Sector Number"; 
1599   binsTrack[1]=180;       xminTrack[1]=0;        xmaxTrack[1]=18; 
1600   //
1601   axisName[2]="R";   axisTitle[2]="r (cm)"; 
1602   binsTrack[2]=53;       xminTrack[2]=85.;        xmaxTrack[2]=245.; 
1603   //
1604   //
1605   axisName[3]="kZ";      axisTitle[3]="dz/dx"; 
1606   binsTrack[3]=36;       xminTrack[3]=-1.8;        xmaxTrack[3]=1.8; 
1607   //
1608   fClusterDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1609   fClusterDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1610   //
1611   //
1612   //
1613   for (Int_t ivar=0;ivar<2;ivar++){
1614     for (Int_t ivar2=0;ivar2<4;ivar2++){
1615       fClusterDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1616       fClusterDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
1617     }
1618   }
1619
1620 }
1621
1622
1623 void AliTPCcalibAlign::MakeResidualHistosTracklet(){
1624   //
1625   // Make tracklet residual histograms
1626   //
1627   Double_t xminTrack[9], xmaxTrack[9];
1628   Int_t    binsTrack[9];
1629   TString  axisName[9],axisTitle[9];
1630   //
1631   // 0 - delta   of interest
1632   // 1 - global  phi in sector number  as float
1633   // 2 - local   x
1634   // 3 - local   ky
1635   // 4 - local   kz
1636   // 5 - sector  1
1637   // 6 - sector  0
1638   // 7 - z position  0
1639
1640   axisName[0]="delta";   axisTitle[0]="#Delta (cm)"; 
1641   binsTrack[0]=60;       xminTrack[0]=-0.6;        xmaxTrack[0]=0.6; 
1642   //
1643   axisName[1]="phi";   axisTitle[1]="#phi"; 
1644   binsTrack[1]=180;       xminTrack[1]=-TMath::Pi();        xmaxTrack[1]=TMath::Pi(); 
1645   //
1646   axisName[2]="localX";   axisTitle[2]="x (cm)"; 
1647   binsTrack[2]=10;       xminTrack[2]=120.;        xmaxTrack[2]=200.; 
1648   //
1649   axisName[3]="kY";      axisTitle[3]="dy/dx"; 
1650   binsTrack[3]=10;       xminTrack[3]=-0.5;        xmaxTrack[3]=0.5; 
1651   //
1652   axisName[4]="kZ";      axisTitle[4]="dz/dx"; 
1653   binsTrack[4]=22;       xminTrack[4]=-1.1;        xmaxTrack[4]=1.1; 
1654   //
1655   axisName[5]="is1";      axisTitle[5]="is1"; 
1656   binsTrack[5]=72;       xminTrack[5]=0;        xmaxTrack[5]=72;
1657   //
1658   axisName[6]="is0";      axisTitle[6]="is0"; 
1659   binsTrack[6]=72;       xminTrack[6]=0;        xmaxTrack[6]=72;
1660   //
1661   axisName[7]="z";        axisTitle[7]="z(cm)"; 
1662   binsTrack[7]=12;        xminTrack[7]=-240;        xmaxTrack[7]=240;
1663   //
1664   axisName[8]="IsPrimary";        axisTitle[8]="Is Primary"; 
1665   binsTrack[8]=2;        xminTrack[8]=-0.1;        xmaxTrack[8]=1.1;
1666  
1667   //
1668   xminTrack[0]=-0.25;        xmaxTrack[0]=0.25; 
1669   fTrackletDelta[0] = new THnSparseF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
1670   xminTrack[0]=-0.5;        xmaxTrack[0]=0.5; 
1671   fTrackletDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
1672   xminTrack[0]=-0.005;        xmaxTrack[0]=0.005; 
1673   fTrackletDelta[2] = new THnSparseF("#Delta_{kY}","#Delta_{kY}", 9, binsTrack,xminTrack, xmaxTrack);
1674   xminTrack[0]=-0.008;        xmaxTrack[0]=0.008; 
1675   fTrackletDelta[3] = new THnSparseF("#Delta_{kZ}","#Delta_{kZ}", 9, binsTrack,xminTrack, xmaxTrack);
1676   //
1677   //
1678   //
1679   for (Int_t ivar=0;ivar<4;ivar++){
1680     for (Int_t ivar2=0;ivar2<9;ivar2++){
1681       fTrackletDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1682       fTrackletDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
1683     }
1684   }
1685
1686 }
1687
1688
1689
1690 void AliTPCcalibAlign::FillHisto(const Double_t *t1,
1691                                  const Double_t *t2,
1692                                  Int_t s1,Int_t s2) {
1693   //
1694   // Fill residual histograms
1695   // Track2-Track1
1696   // Innner-Outer
1697   // Left right - x-y
1698   // A-C side
1699   if (1)  {  
1700     Double_t dy     = t2[1]-t1[1];
1701     Double_t dphi   = t2[2]-t1[2];
1702     Double_t dz     = t2[3]-t1[3];
1703     Double_t dtheta = t2[4]-t1[4];
1704     Double_t zmean = (t2[3]+t1[3])*0.5;
1705     //
1706     GetHisto(kPhi,s1,s2,kTRUE)->Fill(dphi);    
1707     GetHisto(kTheta,s1,s2,kTRUE)->Fill(dtheta);
1708     GetHisto(kY,s1,s2,kTRUE)->Fill(dy);
1709     GetHisto(kZ,s1,s2,kTRUE)->Fill(dz);
1710     //
1711     GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(zmean,dphi);    
1712     GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(zmean,dtheta);
1713     GetHisto(kYz,s1,s2,kTRUE)->Fill(zmean,dy);
1714     GetHisto(kZz,s1,s2,kTRUE)->Fill(zmean,dz);
1715     //
1716     GetHisto(kYPhi,s1,s2,kTRUE)->Fill(t2[2],dy);
1717     GetHisto(kZTheta,s1,s2,kTRUE)->Fill(t2[4],dz);
1718   }     
1719 }
1720
1721
1722 void AliTPCcalibAlign::FillHisto(AliExternalTrackParam *tp1,
1723                                  AliExternalTrackParam *tp2,
1724                                  Int_t s1,Int_t s2) {
1725   //
1726   // Fill residual histograms
1727   // Track2-Track1
1728   if (s2<s1) return;//
1729   const Double_t kEpsilon=0.001;
1730   Double_t x[9]={0,0,0,0,0,0,0,0,0};
1731   AliExternalTrackParam p1(*tp1);
1732   AliExternalTrackParam p2(*tp2);
1733   if (s1%18==s2%18) {
1734     // inner outer - match at the IROC-OROC boundary
1735     if (!p1.PropagateTo(fXIO, AliTrackerBase::GetBz())) return;
1736   }
1737   if (!p2.Rotate(p1.GetAlpha())) return;
1738   if (!p2.PropagateTo(p1.GetX(),AliTrackerBase::GetBz())) return;
1739   if (TMath::Abs(p1.GetX()-p2.GetX())>kEpsilon) return;
1740   Double_t xyz[3];
1741   p1.GetXYZ(xyz);
1742   x[1]=TMath::ATan2(xyz[1],xyz[0]);
1743   x[2]=p1.GetX();
1744   x[3]=0.5*(p1.GetSnp()+p2.GetSnp());  // mean snp
1745   x[4]=0.5*(p1.GetTgl()+p2.GetTgl());  // mean tgl
1746   x[5]=s2;
1747   x[6]=s1;
1748   x[7]=0.5*(p1.GetZ()+p2.GetZ());
1749   // is primary ?
1750   Int_t  isPrimary = (TMath::Abs(p1.GetTgl()-p1.GetZ()/p1.GetX())<0.1) ? 1:0;
1751   x[8]= isPrimary;
1752   //
1753   x[0]=p2.GetY()-p1.GetY();
1754   fTrackletDelta[0]->Fill(x);
1755   x[0]=p2.GetZ()-p1.GetZ();
1756   fTrackletDelta[1]->Fill(x);
1757   x[0]=p2.GetSnp()-p1.GetSnp();
1758   fTrackletDelta[2]->Fill(x);
1759   x[0]=p2.GetTgl()-p1.GetTgl();
1760   fTrackletDelta[3]->Fill(x);
1761   TTreeSRedirector *cstream = GetDebugStreamer();    
1762   if (cstream){
1763     (*cstream)<<"trackletMatch"<<
1764       "tp1.="<<tp1<<   // input tracklet
1765       "tp2.="<<tp2<<
1766       "p1.="<<&p1<<    // tracklet in the ref frame
1767       "p2.="<<&p2<<
1768       "s1="<<s1<<
1769       "s2="<<s2<<
1770       "\n";
1771   }      
1772   
1773 }
1774
1775
1776
1777 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
1778 {
1779   //
1780   // return specified residual histogram - it is only QA
1781   // if force specified the histogram and given histogram is not existing 
1782   //  new histogram is created
1783   //
1784   if (GetIndex(s1,s2)>=72*72) return 0;
1785   TObjArray *histoArray=0;
1786   switch (type) {
1787   case kY:
1788     histoArray = &fDyHistArray; break;
1789   case kZ:
1790     histoArray = &fDzHistArray; break;
1791   case kPhi:
1792     histoArray = &fDphiHistArray; break;
1793   case kTheta:
1794     histoArray = &fDthetaHistArray; break;
1795   case kYPhi:
1796     histoArray = &fDyPhiHistArray; break;
1797   case kZTheta:
1798     histoArray = &fDzThetaHistArray; break;
1799   case kYz:
1800     histoArray = &fDyZHistArray; break;
1801   case kZz:
1802     histoArray = &fDzZHistArray; break;
1803   case kPhiZ:
1804     histoArray = &fDphiZHistArray; break;
1805   case kThetaZ:
1806     histoArray = &fDthetaZHistArray; break;
1807   }
1808   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
1809   if (histo) return histo;
1810   if (force==kFALSE) return 0; 
1811   //
1812   stringstream name;
1813   stringstream title;
1814   switch (type) {    
1815   case kY:
1816     name<<"hist_y_"<<s1<<"_"<<s2;
1817     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1818     histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.5,0.5); // +/- 5 mm
1819     break;
1820   case kZ:
1821     name<<"hist_z_"<<s1<<"_"<<s2;
1822     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1823     histo = new TH1D(name.str().c_str(),title.str().c_str(),100,-0.3,0.3); // +/- 3 mm
1824     break;
1825   case kPhi:
1826     name<<"hist_phi_"<<s1<<"_"<<s2;
1827     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1828     histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.01,0.01); // +/- 10 mrad
1829     break;
1830   case kTheta:
1831     name<<"hist_theta_"<<s1<<"_"<<s2;
1832     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1833     histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.01,0.01); // +/- 10 mrad
1834     break;
1835     //
1836     //
1837   case kYPhi:
1838     name<<"hist_yphi_"<<s1<<"_"<<s2;
1839     title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1840     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,100,-0.5,0.5); // +/- 5 mm
1841     break;
1842   case kZTheta:
1843     name<<"hist_ztheta_"<<s1<<"_"<<s2;
1844     title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1845     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,100,-0.3,0.3); // +/- 3 mm
1846     break;
1847     //
1848     //
1849     //
1850   case kYz:
1851     name<<"hist_yz_"<<s1<<"_"<<s2;
1852     title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1853     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.5,0.5); // +/- 5 mm
1854     break;
1855   case kZz:
1856     name<<"hist_zz_"<<s1<<"_"<<s2;
1857     title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1858     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.3,0.3); // +/- 3 mm
1859     break;
1860   case kPhiZ:
1861     name<<"hist_phiz_"<<s1<<"_"<<s2;
1862     title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1863     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.01,0.01); // +/- 10 mrad
1864     break;
1865   case kThetaZ:
1866     name<<"hist_thetaz_"<<s1<<"_"<<s2;
1867     title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1868     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.01,0.01); // +/- 10 mrad
1869     break;
1870
1871
1872   }
1873   histo->SetDirectory(0);
1874   histoArray->AddAt(histo,GetIndex(s1,s2));
1875   return histo;
1876 }
1877
1878 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec, 
1879                                            Int_t i0, Int_t i1, FitType type) 
1880 {
1881   //
1882   //
1883   //
1884   TMatrixD mat;
1885   //TObjArray *fitArray=0;
1886   Double_t xsec[1000];
1887   Double_t ysec[1000];
1888   Int_t npoints=0;
1889   for (Int_t isec = sec0; isec<=sec1; isec++){
1890     Int_t isec2 = (isec+dsec)%72;    
1891     switch (type) {
1892     case k6:
1893       GetTransformation6(isec,isec2,mat);break;
1894     case k9:
1895       GetTransformation9(isec,isec2,mat);break;
1896     case k12:
1897       GetTransformation12(isec,isec2,mat);break;
1898     }
1899     xsec[npoints]=isec;
1900     ysec[npoints]=mat(i0,i1);
1901     ++npoints;
1902   }
1903   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1904   Char_t name[1000];
1905   snprintf(name,100,"Mat[%d,%d]  Type=%d",i0,i1,type);
1906   gr->SetName(name);
1907   return gr;
1908 }
1909
1910 void  AliTPCcalibAlign::MakeTree(const char *fname, Int_t minPoints){
1911   //
1912   // make tree with alignment cosntant  -
1913   // For  QA visualization
1914   //
1915   /*
1916     TFile fcalib("CalibObjects.root");
1917     TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
1918     AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
1919     align->EvalFitters();
1920     align->MakeTree("alignTree.root");
1921     TFile falignTree("alignTree.root");
1922     TTree * treeAlign = (TTree*)falignTree.Get("Align");
1923    */
1924   TTreeSRedirector cstream(fname);
1925   for (Int_t s1=0;s1<72;++s1)
1926     for (Int_t s2=0;s2<72;++s2){
1927       TMatrixD m6;
1928       TMatrixD m6FX;
1929       TMatrixD m9;
1930       TMatrixD m12;
1931       TVectorD param6Diff;  // align parameters diff 
1932       TVectorD param6s1(6);    // align parameters sector1 
1933       TVectorD param6s2(6);    // align parameters sector2 
1934
1935       //
1936       //
1937       if (fSectorParamA){
1938         TMatrixD * kpar = fSectorParamA;
1939         TMatrixD * kcov = fSectorCovarA;
1940         if (s1%36>=18){
1941           kpar = fSectorParamC;
1942           kcov = fSectorCovarC;
1943         }
1944         for (Int_t ipar=0;ipar<6;ipar++){
1945           Int_t isec1 = s1%18;
1946           Int_t isec2 = s2%18;
1947           if (s1>35) isec1+=18;
1948           if (s2>35) isec2+=18; 
1949           param6s1(ipar)=(*kpar)(6*isec1+ipar,0);
1950           param6s2(ipar)=(*kpar)(6*isec2+ipar,0);
1951         }
1952       }
1953         
1954       Double_t dy=0, dz=0, dphi=0,dtheta=0;
1955       Double_t sy=0, sz=0, sphi=0,stheta=0;
1956       Double_t ny=0, nz=0, nphi=0,ntheta=0;
1957       Double_t chi2v12=0, chi2v9=0, chi2v6=0;
1958       //      Int_t npoints=0;
1959       // TLinearFitter * fitter = 0;      
1960       if (fPoints[GetIndex(s1,s2)]>minPoints){
1961         //
1962         //
1963         //
1964 //      fitter = GetFitter12(s1,s2);
1965 //      npoints = fitter->GetNpoints();
1966 //      chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1967         
1968 //      //
1969 //      fitter = GetFitter9(s1,s2);
1970 //      npoints = fitter->GetNpoints();
1971 //      chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1972 //      //
1973 //      fitter = GetFitter6(s1,s2);
1974 //      npoints = fitter->GetNpoints();
1975 //      chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1976 //      fitter->GetParameters(param6Diff);
1977 //      //
1978 //      GetTransformation6(s1,s2,m6);
1979 //      GetTransformation9(s1,s2,m9);
1980 //      GetTransformation12(s1,s2,m12);
1981 //      //
1982 //      fitter = GetFitter6(s1,s2);
1983 //      //fitter->FixParameter(3,0);
1984 //      //fitter->Eval();
1985 //      GetTransformation6(s1,s2,m6FX);
1986         //
1987         TH1 * his=0;
1988         his = GetHisto(kY,s1,s2);
1989         if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1990         his = GetHisto(kZ,s1,s2);
1991         if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1992         his = GetHisto(kPhi,s1,s2);
1993         if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
1994         his = GetHisto(kTheta,s1,s2);
1995         if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
1996         //
1997
1998       }
1999       AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2000       if (!magF) AliError("Magneticd field - not initialized");
2001       Double_t bz = magF->SolenoidField()/10.; //field in T
2002
2003       cstream<<"Align"<<
2004         "run="<<fRun<<  // run
2005         "bz="<<bz<<
2006         "s1="<<s1<<     // reference sector
2007         "s2="<<s2<<     // sector to align
2008         "m6FX.="<<&m6FX<<   // tranformation matrix
2009         "m6.="<<&m6<<   // tranformation matrix
2010         "m9.="<<&m9<<   // 
2011         "m12.="<<&m12<<
2012         "chi2v12="<<chi2v12<<
2013         "chi2v9="<<chi2v9<<
2014         "chi2v6="<<chi2v6<<
2015         //
2016         "p6.="<<&param6Diff<<
2017         "p6s1.="<<&param6s1<<
2018         "p6s2.="<<&param6s2<<
2019         //               histograms mean RMS and entries
2020         "dy="<<dy<<  
2021         "sy="<<sy<<
2022         "ny="<<ny<<
2023         "dz="<<dz<<
2024         "sz="<<sz<<
2025         "nz="<<nz<<
2026         "dphi="<<dphi<<
2027         "sphi="<<sphi<<
2028         "nphi="<<nphi<<
2029         "dtheta="<<dtheta<<
2030         "stheta="<<stheta<<
2031         "ntheta="<<ntheta<<
2032         "\n";
2033     }
2034
2035 }
2036
2037
2038 //_____________________________________________________________________
2039 Long64_t AliTPCcalibAlign::Merge(TCollection* const list) {
2040   //
2041   // merge function 
2042   //
2043   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
2044   if (!list)
2045     return 0;  
2046   if (list->IsEmpty())
2047     return 1;
2048   
2049   TIterator* iter = list->MakeIterator();
2050   TObject* obj = 0;  
2051   iter->Reset();
2052   Int_t count=0;
2053   //
2054   TString str1(GetName());
2055   while((obj = iter->Next()) != 0)
2056     {
2057       AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
2058       if (entry == 0) continue; 
2059       if (str1.CompareTo(entry->GetName())!=0) continue;
2060       Add(entry);
2061       count++;
2062     } 
2063   return count;
2064 }
2065
2066
2067 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
2068   //
2069   // Add entry - used for merging of compoents
2070   //
2071   for (Int_t i=0; i<72;i++){
2072     for (Int_t j=0; j<72;j++){
2073       if (align->fPoints[GetIndex(i,j)]<1) continue;
2074       fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
2075       //
2076       //
2077       //
2078       for (Int_t itype=0; itype<10; itype++){
2079         TH1 * his0=0, *his1=0;
2080         his0 = GetHisto((HistoType)itype,i,j);
2081         his1 = align->GetHisto((HistoType)itype,i,j);
2082         if (his1){
2083           if (his0) his0->Add(his1);
2084           else {
2085             his0 = GetHisto((HistoType)itype,i,j,kTRUE);
2086             his0->Add(his1);
2087           }
2088         }       
2089       }      
2090     }
2091   }
2092   TLinearFitter *f0=0;
2093   TLinearFitter *f1=0;
2094   for (Int_t i=0; i<72;i++){
2095     for (Int_t j=0; j<72;j++){     
2096       if (align->fPoints[GetIndex(i,j)]<1) continue;
2097       //
2098       //
2099       // fitter12
2100       f0 =  GetFitter12(i,j);
2101       f1 =  align->GetFitter12(i,j);
2102       if (f1){
2103         if (f0) f0->Add(f1);
2104         else {
2105           fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
2106         }
2107       }      
2108       //
2109       // fitter9
2110       f0 =  GetFitter9(i,j);
2111       f1 =  align->GetFitter9(i,j);
2112       if (f1){
2113         if (f0) f0->Add(f1);
2114         else { 
2115           fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
2116         }
2117       }      
2118       f0 =  GetFitter6(i,j);
2119       f1 =  align->GetFitter6(i,j);
2120       if (f1){
2121         if (f0) f0->Add(f1);
2122         else {
2123           fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
2124         }
2125       }   
2126     }
2127   }
2128   //
2129   // Add Kalman filter
2130   //
2131   for (Int_t i=0;i<36;i++){
2132     TMatrixD *par0 = (TMatrixD*)fArraySectorIntParam.At(i);
2133     if (!par0){
2134       MakeSectorKalman();
2135       par0 = (TMatrixD*)fArraySectorIntParam.At(i);      
2136     }
2137     TMatrixD *par1 = (TMatrixD*)align->fArraySectorIntParam.At(i);
2138     if (!par1) continue;
2139     //
2140     TMatrixD *cov0 = (TMatrixD*)fArraySectorIntCovar.At(i);
2141     TMatrixD *cov1 = (TMatrixD*)align->fArraySectorIntCovar.At(i);
2142     UpdateSectorKalman(*par0,*cov0,*par1,*cov1);
2143   }
2144   if (!fSectorParamA){
2145     MakeKalman();
2146   }
2147   if (align->fSectorParamA){
2148     UpdateKalman(*fSectorParamA,*fSectorCovarA,*align->fSectorParamA,*align->fSectorCovarA);
2149     UpdateKalman(*fSectorParamC,*fSectorCovarC,*align->fSectorParamC,*align->fSectorCovarC);
2150   }
2151   if (!fClusterDelta[0]) MakeResidualHistos();
2152
2153   for (Int_t i=0; i<2; i++){
2154     if (align->fClusterDelta[i]){
2155       fClusterDelta[i]->Add(align->fClusterDelta[i]);
2156   //     align->fClusterDelta[i]->GetAxis(0)->SetRangeUser(-0.87,0.87);
2157 //       align->fClusterDelta[i]->GetAxis(3)->SetRangeUser(-0.87,0.87);
2158 //       fClusterDelta[i]->GetAxis(0)->SetRangeUser(-0.87,0.87);
2159 //       fClusterDelta[i]->GetAxis(3)->SetRangeUser(-0.87,0.87);
2160 //       Int_t idim[4]={0,1,2,3};
2161 //       THnSparse *htemp=align->fClusterDelta[i]->Projection(4,idim);
2162 //       THnSparse *htemp1=fClusterDelta[i]->Projection(4,idim);      
2163 //       htemp1->Add(htemp);
2164 //       delete fClusterDelta[i];
2165 //       fClusterDelta[i]=htemp1;      
2166 //       delete htemp;
2167     }
2168   }
2169   
2170   for (Int_t i=0; i<4; i++){
2171     if (!fTrackletDelta[i] && align->fTrackletDelta[i]) {
2172       fTrackletDelta[i]= (THnSparse*)(align->fTrackletDelta[i]->Clone());
2173       continue;
2174     }
2175     if (align->fTrackletDelta[i]) {
2176       fTrackletDelta[i]->Add(align->fTrackletDelta[i]);
2177       //
2178    //    align->fTrackletDelta[i]->GetAxis(3)->SetRangeUser(-0.36,0.36);
2179 //       align->fTrackletDelta[i]->GetAxis(4)->SetRangeUser(-0.87,0.87);
2180 //       fTrackletDelta[i]->GetAxis(3)->SetRangeUser(-0.36,0.36);
2181 //       fTrackletDelta[i]->GetAxis(4)->SetRangeUser(-0.87,0.87);
2182 //       //
2183 //       Int_t idim[9]={0,1,2,3,4,5,6,7,8};
2184 //       THnSparse *htemp=align->fTrackletDelta[i]->Projection(9,idim);
2185 //       THnSparse *htemp1=fTrackletDelta[i]->Projection(9,idim);
2186 //       htemp1->Add(htemp);
2187 //       delete fTrackletDelta[i];
2188 //       fTrackletDelta[i]=htemp1;
2189 //       delete htemp;
2190     }
2191   }
2192
2193 }
2194
2195 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){
2196   //
2197   // GetTransformed value
2198   //
2199   //
2200   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
2201   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
2202   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
2203   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2204   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2205   
2206   
2207   const TMatrixD * mat = GetTransformation(s1,s2,type);
2208   if (!mat) {
2209     if (value==0) return x1;
2210     if (value==1) return y1;
2211     if (value==2) return z1;
2212     if (value==3) return dydx1;
2213     if (value==4) return dzdx1;
2214     //
2215     if (value==5) return dydx1;
2216     if (value==6) return dzdx1;
2217     return 0;
2218   }
2219   Double_t valT=0;
2220
2221   if (value==0){
2222     valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
2223   }
2224
2225   if (value==1){
2226     valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
2227   }
2228   if (value==2){
2229     valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
2230   }
2231   if (value==3){
2232     //    dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2233     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
2234     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
2235   }
2236
2237   if (value==4){
2238     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)    
2239     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
2240     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
2241   }
2242   //
2243   if (value==5){
2244     // onlys shift in angle
2245     //    dydx2 =  (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
2246     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1;
2247   }
2248
2249   if (value==6){
2250     // only shift in angle
2251     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)    
2252     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1;
2253   }
2254   //
2255   return valT;
2256 }
2257
2258
2259 void  AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){
2260   //
2261   // Update track parameters t1
2262   //
2263   TMatrixD vecXk(5,1);    // X vector
2264   TMatrixD covXk(5,5);    // X covariance 
2265   TMatrixD matHk(1,5);    // vector to mesurement
2266   TMatrixD measR(1,1);    // measurement error 
2267   //TMatrixD matQk(5,5);    // prediction noise vector
2268   TMatrixD vecZk(1,1);    // measurement
2269   //
2270   TMatrixD vecYk(1,1);    // Innovation or measurement residual
2271   TMatrixD matHkT(5,1);
2272   TMatrixD matSk(1,1);    // Innovation (or residual) covariance
2273   TMatrixD matKk(5,1);    // Optimal Kalman gain
2274   TMatrixD mat1(5,5);     // update covariance matrix
2275   TMatrixD covXk2(5,5);   // 
2276   TMatrixD covOut(5,5);
2277   //
2278   Double_t *param1=(Double_t*) track1.GetParameter();
2279   Double_t *covar1=(Double_t*) track1.GetCovariance();
2280
2281   //
2282   // copy data to the matrix
2283   for (Int_t ipar=0; ipar<5; ipar++){
2284     vecXk(ipar,0) = param1[ipar];
2285     for (Int_t jpar=0; jpar<5; jpar++){
2286       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
2287     }
2288   }
2289   //
2290   //
2291   //
2292   vecZk(0,0) = track2.GetParameter()[4];   // 1/pt measurement from track 2
2293   measR(0,0) = track2.GetCovariance()[14];  // 1/pt measurement error
2294   if (noField) {
2295     measR(0,0)*=0.000000001;
2296     vecZk(0,0)=0.;
2297   }
2298   //
2299   matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0;  
2300   matHk(0,3)= 0;    matHk(0,4)= 1;           // vector to measurement
2301   //
2302   //
2303   //
2304   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
2305   matHkT=matHk.T(); matHk.T();
2306   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
2307   matSk.Invert();
2308   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
2309   vecXk += matKk*vecYk;                      //  updated vector 
2310   mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
2311   covXk2 = (mat1-(matKk*matHk));
2312   covOut =  covXk2*covXk; 
2313   //
2314   //
2315   //
2316   // copy from matrix to parameters
2317   if (0) {
2318     covOut.Print();
2319     vecXk.Print();
2320     covXk.Print();
2321     track1.Print();
2322     track2.Print();
2323   }
2324
2325   for (Int_t ipar=0; ipar<5; ipar++){
2326     param1[ipar]= vecXk(ipar,0) ;
2327     for (Int_t jpar=0; jpar<5; jpar++){
2328       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
2329     }
2330   }
2331
2332 }
2333
2334 void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){
2335   //
2336   //  Global Align -combine the partial alignment of pair of sectors
2337   //  minPoints - minimal number of points - don't use sector alignment wit smaller number
2338   //  sysError  - error added to the alignemnt error
2339   //
2340   AliTPCcalibAlign * align = this;
2341   TMatrixD * arrayAlign[72];
2342   TMatrixD * arrayAlignDiff[72];
2343   //
2344   for (Int_t i=0;i<72; i++) {
2345     TMatrixD * mat = new TMatrixD(4,4);
2346     mat->UnitMatrix();
2347     arrayAlign[i]=mat;
2348     arrayAlignDiff[i]=(TMatrixD*)(mat->Clone());
2349   }
2350
2351   TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root");
2352   for (Int_t iter=0; iter<niter;iter++){
2353     printf("Iter=\t%d\n",iter);
2354     for (Int_t is0=0;is0<72; is0++) {
2355       //
2356       //TMatrixD  *mati0 = arrayAlign[is0];
2357       TMatrixD matDiff(4,4);
2358       Double_t sumw=0;      
2359       for (Int_t is1=0;is1<72; is1++) {
2360         Bool_t invers=kFALSE;
2361         Int_t npoints=0;
2362         TMatrixD covar;
2363         TVectorD errors;
2364         const TMatrixD *mat = align->GetTransformation(is0,is1,0); 
2365         if (mat){
2366           npoints = align->GetFitter6(is0,is1)->GetNpoints();
2367           if (npoints>minPoints){
2368             align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar);
2369             align->GetFitter6(is0,is1)->GetErrors(errors);
2370           }
2371         }
2372         else{
2373           invers=kTRUE;
2374           mat = align->GetTransformation(is1,is0,0); 
2375           if (mat) {
2376             npoints = align->GetFitter6(is1,is0)->GetNpoints();
2377             if (npoints>minPoints){
2378               align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar);
2379               align->GetFitter6(is1,is0)->GetErrors(errors);
2380             }
2381           }
2382         }
2383         if (!mat) continue;
2384         if (npoints<minPoints) continue;
2385         //
2386         Double_t weight=1;
2387         if (is1/36>is0/36) weight*=2./3.; //IROC-OROC
2388         if (is1/36<is0/36) weight*=1./3.; //OROC-IROC
2389         if (is1/36==is0/36) weight*=1/3.; //OROC-OROC
2390         if (is1%36!=is0%36) weight*=1/2.; //Not up-down
2391         weight/=(errors[4]*errors[4]+sysError*sysError);  // wieghting with error in Y
2392         //
2393         //
2394         TMatrixD matT = *mat;   
2395         if (invers)  matT.Invert();
2396         TMatrixD diffMat= (*(arrayAlign[is1]))*matT;
2397         diffMat-=(*arrayAlign[is0]);
2398         matDiff+=weight*diffMat;
2399         sumw+=weight;
2400
2401         (*cstream)<<"LAlign"<<
2402           "iter="<<iter<<
2403           "s0="<<is0<<
2404           "s1="<<is1<<
2405           "npoints="<<npoints<<
2406           "m60.="<<arrayAlign[is0]<<
2407           "m61.="<<arrayAlign[is1]<<
2408           "m01.="<<&matT<<
2409           "diff.="<<&diffMat<<
2410           "cov.="<<&covar<<
2411           "err.="<<&errors<<
2412           "\n";
2413       }
2414       if (sumw>0){
2415         matDiff*=1/sumw;
2416         matDiff(0,0)=0;
2417         matDiff(1,1)=0;
2418         matDiff(1,1)=0;
2419         matDiff(1,1)=0; 
2420         (*arrayAlignDiff[is0]) = matDiff;       
2421       }
2422     }
2423     for (Int_t is0=0;is0<72; is0++) {
2424       if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]);
2425       if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]);
2426        //
2427       (*cstream)<<"GAlign"<<
2428         "iter="<<iter<<
2429         "s0="<<is0<<
2430         "m6.="<<arrayAlign[is0]<<
2431         "\n";
2432     }    
2433   }
2434
2435   delete cstream;
2436   for (Int_t isec=0;isec<72;isec++){
2437     fCombinedMatrixArray6.AddAt(arrayAlign[isec],isec);
2438     delete arrayAlignDiff[isec];
2439   }
2440 }
2441
2442
2443  Int_t  AliTPCcalibAlign::RefitLinear(const AliTPCseed * track, Int_t isec, Double_t *fitParam, Int_t refSector,  TMatrixD &tparam, TMatrixD&tcovar, Double_t xRef, Bool_t both){
2444   //
2445   // Refit tracklet linearly using clusters at  given sector isec
2446   // Clusters are rotated to the  reference frame of sector refSector
2447   // 
2448   // fit parameters and errors retruning in the fitParam
2449   //
2450   // seed      - acces to the original clusters
2451   // isec      - sector to be refited
2452   // fitParam  - 
2453   //           0  lx             
2454   //           1  ly
2455   //           2  dy/dz
2456   //           3  lz
2457   //           4  dz/dx
2458   //           5  sx 
2459   //           6  sy
2460   //           7  sdydx
2461   //           8  sz
2462   //           9  sdzdx
2463   // ref sector is the sector defining ref frame - rotation
2464   // return value - number of used clusters
2465
2466   const Int_t      kMinClusterF=15;
2467   const  Int_t     kdrow1 =10;        // rows to skip at the end      
2468   const  Int_t     kdrow0 =3;        // rows to skip at beginning  
2469   const  Float_t   kedgeyIn=2.5;
2470   const  Float_t   kedgeyOut=4.0;
2471   const  Float_t   kMaxDist=5;       // max distance -in sigma 
2472   const  Float_t   kMaxCorrY=0.05;    // max correction
2473   //
2474   Double_t dalpha = 0;
2475   if ((refSector%18)!=(isec%18)){
2476     dalpha = -((refSector%18)-(isec%18))*TMath::TwoPi()/18.;
2477   }
2478   Double_t ca = TMath::Cos(dalpha);
2479   Double_t sa = TMath::Sin(dalpha);
2480   //
2481   //
2482   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
2483   //
2484   // full track fit parameters
2485   // 
2486   static TLinearFitter fyf(2,"pol1");   // change to static - suggestion of calgrind - 30 % of time
2487   static TLinearFitter fzf(2,"pol1");   // relative to time of given class
2488   TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
2489   TMatrixD covY(4,4),covZ(4,4);
2490   Double_t chi2FacY =1;
2491   Double_t chi2FacZ =1;
2492   Int_t nf=0;
2493   //
2494   //
2495   //
2496   Float_t erry=0.1;   // initial cluster error estimate
2497   Float_t errz=0.1;   // initial cluster error estimate
2498   for (Int_t iter=0; iter<2; iter++){
2499     fyf.ClearPoints();
2500     fzf.ClearPoints();
2501     for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
2502       AliTPCclusterMI *c=track->GetClusterPointer(irow);
2503       if (!c) continue;      
2504       //
2505       if (c->GetDetector()%36!=(isec%36)) continue;
2506       if (!both && c->GetDetector()!=isec) continue;
2507
2508       if (c->GetRow()<kdrow0) continue;
2509       //cluster position in reference frame 
2510       Double_t lxR   =   ca*c->GetX()-sa*c->GetY();  
2511       Double_t lyR   =  +sa*c->GetX()+ca*c->GetY();
2512       Double_t lzR   =  c->GetZ();
2513
2514       Double_t dx = lxR -xRef;   // distance to reference X
2515       Double_t x[2]={dx, dx*dx};
2516
2517       Double_t yfitR  =    pyf[0]+pyf[1]*dx;  // fit value Y in ref frame
2518       Double_t zfitR  =    pzf[0]+pzf[1]*dx;  // fit value Z in ref frame
2519       //
2520       Double_t yfit   =  -sa*lxR + ca*yfitR;  // fit value Y in local frame
2521       //
2522       if (iter==0 &&c->GetType()<0) continue;
2523       if (iter>0){      
2524         if (TMath::Abs(lyR-yfitR)>kMaxDist*erry) continue;
2525         if (TMath::Abs(lzR-zfitR)>kMaxDist*errz) continue;
2526         Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2527         if (isec<36 && dedge<kedgeyIn) continue;
2528         if (isec>35 && dedge<kedgeyOut) continue;
2529         Double_t corrtrY =  
2530           corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
2531                                   c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2532         Double_t corrclY =  
2533           corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
2534                                   c->GetY(),c->GetY(), c->GetZ(), pyf[1], c->GetMax(),2.5);
2535         if (TMath::Abs((corrtrY+corrclY)*0.5)>kMaxCorrY) continue;
2536         if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2537       }
2538       fyf.AddPoint(x,lyR,erry);
2539       fzf.AddPoint(x,lzR,errz);
2540     }
2541     nf = fyf.GetNpoints();
2542     if (nf<kMinClusterF) return 0;   // not enough points - skip 
2543     fyf.Eval(); 
2544     fyf.GetParameters(pyf); 
2545     fyf.GetErrors(peyf);
2546     fzf.Eval(); 
2547     fzf.GetParameters(pzf); 
2548     fzf.GetErrors(pezf);
2549     chi2FacY = TMath::Sqrt(fyf.GetChisquare()/(fyf.GetNpoints()-2.));
2550     chi2FacZ = TMath::Sqrt(fzf.GetChisquare()/(fzf.GetNpoints()-2.));
2551     peyf[0]*=chi2FacY;
2552     peyf[1]*=chi2FacY;
2553     pezf[0]*=chi2FacZ;
2554     pezf[1]*=chi2FacZ;
2555     erry*=chi2FacY;
2556     errz*=chi2FacZ;
2557     fyf.GetCovarianceMatrix(covY);
2558     fzf.GetCovarianceMatrix(covZ);
2559     for (Int_t i0=0;i0<2;i0++)
2560       for (Int_t i1=0;i1<2;i1++){
2561         covY(i0,i1)*=chi2FacY*chi2FacY;
2562         covZ(i0,i1)*=chi2FacZ*chi2FacZ;
2563       }
2564   }
2565   fitParam[0] = xRef;
2566   //
2567   fitParam[1] = pyf[0];
2568   fitParam[2] = pyf[1];
2569   fitParam[3] = pzf[0];
2570   fitParam[4] = pzf[1];
2571   //
2572   fitParam[5] = 0;
2573   fitParam[6] = peyf[0];
2574   fitParam[7] = peyf[1];
2575   fitParam[8] = pezf[0];
2576   fitParam[9] = pezf[1];
2577   //
2578   //
2579   tparam(0,0) = pyf[0];
2580   tparam(1,0) = pyf[1];
2581   tparam(2,0) = pzf[0];
2582   tparam(3,0) = pzf[1];
2583   //
2584   tcovar(0,0) = covY(0,0);
2585   tcovar(1,1) = covY(1,1);
2586   tcovar(1,0) = covY(1,0);
2587   tcovar(0,1) = covY(0,1); 
2588   tcovar(2,2) = covZ(0,0);
2589   tcovar(3,3) = covZ(1,1);
2590   tcovar(3,2) = covZ(1,0);
2591   tcovar(2,3) = covZ(0,1);
2592   return nf;
2593 }
2594
2595 void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){
2596   //
2597   // Update the cluster residula histograms for setup with field
2598   // Kalman track fitting is used
2599   //  Only high momenta primary tracks used
2600   //
2601   // 1. Apply selection
2602   // 2. Refit the track - in-out
2603   // 3. Refit the track - out-in
2604   // 4. Combine In and Out track - - fil cluster residuals
2605   //
2606   const Double_t kPtCut=1.0;    // pt
2607   const Double_t kSnpCut=0.2; // snp cut
2608   const Double_t kNclCut=120; //
2609   const Double_t kVertexCut=1;
2610   const Double_t kMaxDist=0.5; // max distance between tracks and cluster
2611   const Double_t kEdgeCut = 2.5;
2612   const Double_t kDelta2=0.2*0.2;  // initial increase in covar matrix
2613   const Double_t kSigma=0.3;       // error increase towards edges of TPC 
2614   const Double_t kSkipBoundary=7.5;  // skip track updates in the boundary IFC,OFC, IO
2615   //
2616   if (!fCurrentTrack) return;
2617   if (!fCurrentFriendTrack) return;
2618   Float_t vertexXY=0,vertexZ=0;
2619   fCurrentTrack->GetImpactParameters(vertexXY,vertexZ);
2620   if (TMath::Abs(vertexXY)>kVertexCut) return;
2621   if (TMath::Abs(vertexZ)>kVertexCut) return;
2622   if (TMath::Abs(seed->Pt())<kPtCut) return;
2623   if (seed->GetNumberOfClusters()<kNclCut) return;
2624   if (TMath::Abs(seed->GetSnp())>kSnpCut) return; 
2625   if (!fClusterDelta[0])  MakeResidualHistos();
2626   //
2627   AliExternalTrackParam fitIn[160];
2628   AliExternalTrackParam fitOut[160];
2629   AliTPCROC * roc = AliTPCROC::Instance();
2630   Double_t xmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
2631   Double_t xDiff     = ( -roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
2632   Double_t xIFC      = ( roc->GetPadRowRadii(0,0));
2633   Double_t xOFC      = ( roc->GetPadRowRadii(36,roc->GetNRows(36)-1));
2634   //
2635   Int_t detector=-1;
2636   //
2637   //
2638   AliExternalTrackParam trackIn  = *(fCurrentTrack->GetInnerParam());
2639   AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut());
2640   trackIn.ResetCovariance(10);
2641   trackOut.ResetCovariance(10);
2642   Double_t *covarIn = (Double_t*)trackIn.GetCovariance();
2643   Double_t *covarOut = (Double_t*)trackOut.GetCovariance();
2644   covarIn[0]+=kDelta2; covarIn[2]+=kDelta2; 
2645   covarIn[5]+=kDelta2/(100.*100.); covarIn[9]=kDelta2/(100.*100.);
2646   covarIn[14]+=kDelta2/(5.*5.);
2647   covarOut[0]+=kDelta2; covarOut[2]+=kDelta2; 
2648   covarOut[5]+=kDelta2/(100.*100.); covarOut[9]=kDelta2/(100.*100.);
2649   covarOut[14]+=kDelta2/(5.*5.);
2650   //
2651   static Double_t mass =    TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2652   //  
2653   Int_t ncl=0;
2654   for (Int_t irow=0; irow<160; irow++){
2655     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2656     if (!cl) continue;
2657     if (cl->GetX()<80) continue;
2658     if (detector<0) detector=cl->GetDetector()%36;
2659     if (detector!=cl->GetDetector()%36) return; // cluster from different sectors
2660     // skip such tracks
2661     ncl++;
2662   }
2663   if (ncl<kNclCut) return;
2664   Int_t nclIn=0,nclOut=0;
2665   Double_t xyz[3];
2666   //
2667   // Refit out - store residual maps
2668   //
2669   for (Int_t irow=0; irow<160; irow++){
2670     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2671     if (!cl) continue;
2672     if (cl->GetX()<80) continue;
2673     if (detector<0) detector=cl->GetDetector()%36;
2674     Int_t sector = cl->GetDetector();
2675     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();    
2676     if (cl->GetDetector()%36!=detector) continue;
2677     if (TMath::Abs(dalpha)>0.01){
2678       if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
2679     }
2680     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};    
2681     Double_t cov[3]={0.1,0.,0.1}; 
2682     Double_t dedge =  cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackOut.GetY());
2683     Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
2684     dmiddle*=dmiddle;
2685     //
2686     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2687     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2688     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2689     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2690     cov[0]+=kSigma/dedge;      // bigger error close to the boundary
2691     cov[2]+=kSigma/dedge;      // bigger error close to the boundary
2692     cov[0]*=cov[0];
2693     cov[2]*=cov[2];
2694     if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;    
2695     if (TMath::Abs(dedge)<kEdgeCut) continue;
2696     //
2697     Bool_t doUpdate=kTRUE;
2698     if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
2699     if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
2700     if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
2701     //
2702     if (TMath::Abs(cl->GetY()-trackOut.GetY())<kMaxDist){
2703       nclOut++;
2704       if (doUpdate) trackOut.Update(&r[1],cov);
2705     }
2706     fitOut[irow]=trackOut;
2707   }
2708   
2709   //
2710   // Refit In - store residual maps
2711   //
2712   for (Int_t irow=159; irow>=0; irow--){
2713     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2714     if (!cl) continue;
2715     if (cl->GetX()<80) continue;
2716     if (detector<0) detector=cl->GetDetector()%36;
2717     Int_t sector = cl->GetDetector();
2718     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();    
2719     if (cl->GetDetector()%36!=detector) continue;
2720     if (TMath::Abs(dalpha)>0.01){
2721       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
2722     }
2723     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};    
2724     Double_t cov[3]={0.1,0.,0.1}; 
2725     Double_t dedge =  cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackIn.GetY());
2726     Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
2727     dmiddle*=dmiddle;
2728     //
2729     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2730     cov[0]+=kSigma*dmiddle;     // bigger error at boundary
2731     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2732     cov[2]+=kSigma*dmiddle;     // bigger error at boundary
2733     cov[0]+=kSigma/dedge;      // bigger error close to the boundary
2734     cov[2]+=kSigma/dedge;      // bigger error close to the boundary
2735     cov[0]*=cov[0];
2736     cov[2]*=cov[2];
2737     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;    
2738     if (TMath::Abs(dedge)<kEdgeCut) continue;
2739     Bool_t doUpdate=kTRUE;
2740     if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
2741     if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
2742     if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
2743     if (TMath::Abs(cl->GetY()-trackIn.GetY())<kMaxDist){
2744       nclIn++;
2745       if (doUpdate) trackIn.Update(&r[1],cov);
2746     }
2747     fitIn[irow]=trackIn;
2748   }
2749   //
2750   //
2751   for (Int_t irow=159; irow>=0; irow--){
2752     //
2753     // Update kalman - +- direction
2754     // Store cluster residuals
2755     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2756     if (!cl) continue;
2757     if (cl->GetX()<80) continue;
2758     if (detector<0) detector=cl->GetDetector()%36;
2759     if (cl->GetDetector()%36!=detector) continue;
2760     AliExternalTrackParam trackSmooth = fitIn[irow];
2761     AliTrackerBase::UpdateTrack(trackSmooth, fitOut[irow]);
2762     //
2763     Double_t resVector[5]; 
2764     trackSmooth.GetXYZ(xyz);
2765     resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
2766     if (resVector[1]<0) resVector[1]+=18;
2767     resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
2768     resVector[3]= cl->GetZ()/resVector[2];
2769     //
2770     resVector[0]= cl->GetY()-trackSmooth.GetY();
2771     fClusterDelta[0]->Fill(resVector);
2772     resVector[0]= cl->GetZ()-trackSmooth.GetZ();
2773     fClusterDelta[1]->Fill(resVector);
2774   }
2775
2776 }
2777
2778
2779 void  AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
2780   //
2781   // Update Kalman filter of Alignment - only setup without filed
2782   //       IROC - OROC quadrants
2783   //
2784   if (TMath::Abs(AliTracker::GetBz())>0.5) return;
2785   if (!fClusterDelta[0])  MakeResidualHistos();
2786   //  const Int_t kMinClusterF=40;
2787   const Int_t kMinClusterFit=10;
2788   const Int_t kMinClusterQ=10;
2789   //
2790   const  Int_t     kdrow1Fit =5;         // rows to skip from fit at the end      
2791   const  Int_t     kdrow0Fit =10;        // rows to skip from fit at beginning  
2792   const  Float_t   kedgey=3.0;
2793   const  Float_t   kMaxDist=1;
2794   const  Float_t   kMaxCorrY=0.05;
2795   const  Float_t   kPRFWidth = 0.6;   //cut 2 sigma of PRF
2796   isec = isec%36;     // use the hardware numbering
2797   //
2798   //
2799   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
2800   //
2801   // full track fit parameters
2802   // 
2803   static TLinearFitter fyf(2,"pol1");   // make it static - too much time for comiling of formula
2804   static TLinearFitter fzf(2,"pol1");   // calgrind recomendation
2805   TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
2806   TVectorD pyfc(2),pzfc(2);
2807   Int_t nf=0;
2808   //
2809   // make full fit as reference
2810   //
2811   for (Int_t iter=0; iter<2; iter++){
2812     fyf.ClearPoints();
2813     fzf.ClearPoints();
2814     for (Int_t irow=kdrow0Fit;irow<159-kdrow1Fit;irow++) {
2815       AliTPCclusterMI *c=track->GetClusterPointer(irow);
2816       if (!c) continue;
2817       if ((c->GetDetector()%36)!=isec) continue;
2818       if (c->GetRow()<kdrow0Fit) continue;
2819       Double_t dx = c->GetX()-fXmiddle;
2820       Double_t x[2]={dx, dx*dx};
2821       if (iter==0 &&c->GetType()<0) continue;
2822       if (iter==1){     
2823         Double_t yfit  =  pyf[0]+pyf[1]*dx;
2824         Double_t zfit  =  pzf[0]+pzf[1]*dx;
2825         Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2826         if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
2827         if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
2828         if (dedge<kedgey) continue;
2829         Double_t corrtrY =  
2830           corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
2831                                   c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2832         if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2833       }
2834       if (TMath::Abs(x[0])<10){
2835         fyf.AddPoint(x,c->GetY(),0.1); //use only middle rows+-10cm
2836         fzf.AddPoint(x,c->GetZ(),0.1);            
2837       }
2838     }
2839     nf = fyf.GetNpoints();
2840     if (fyf.GetNpoints()<kMinClusterFit) return;   // not enough points - skip 
2841     if (fzf.GetNpoints()<kMinClusterFit) return;   // not enough points - skip 
2842     fyf.Eval(); 
2843     fyf.GetParameters(pyf); 
2844     fyf.GetErrors(peyf);
2845     fzf.Eval(); 
2846     fzf.GetParameters(pzf); 
2847     fzf.GetErrors(pezf);
2848   }
2849   //
2850   //
2851   //
2852   TVectorD vecX(160);          // x         vector
2853   TVectorD vecY(160);          // residuals vector
2854   TVectorD vecZ(160);                              // residuals vector
2855   TVectorD vPosG(3);                  //vertex position
2856   TVectorD vPosL(3);                 // vertex position in the TPC local system
2857   TVectorF vImpact(2);               //track impact parameter
2858   //  Double_t tofSignal= fCurrentTrack->GetTOFsignal();      // tof signal
2859   TVectorF tpcPosG(3);                                    // global position of track at the middle of fXmiddle
2860   Double_t lphi = TMath::ATan2(pyf[0],fXmiddle);          // expected local phi angle - if vertex at 0
2861   Double_t gphi = 2.*TMath::Pi()*(isec%18+0.5)/18.+lphi;  // expected global phi if vertex at 0
2862   Double_t ky   = pyf[0]/fXmiddle;
2863   Double_t kyE  =0, kzE=0;    // ky and kz expected
2864   Double_t alpha =2.*TMath::Pi()*(isec%18+0.5)/18.;
2865   Double_t scos=TMath::Cos(alpha);
2866   Double_t ssin=TMath::Sin(alpha);
2867   const AliESDVertex* vertex = fCurrentEvent->GetPrimaryVertexTracks();
2868   vertex->GetXYZ(vPosG.GetMatrixArray());
2869   fCurrentTrack->GetImpactParameters(vImpact[0],vImpact[1]);  // track impact parameters
2870   //
2871   tpcPosG[0]= scos*fXmiddle-ssin*pyf[0];
2872   tpcPosG[1]=+ssin*fXmiddle+scos*pyf[0];
2873   tpcPosG[2]=pzf[0];
2874   vPosL[0]= scos*vPosG[0]+ssin*vPosG[1];
2875   vPosL[1]=-ssin*vPosG[0]+scos*vPosG[1];
2876   vPosL[2]= vPosG[2];
2877   kyE = (pyf[0]-vPosL[1])/(fXmiddle-vPosL[0]);
2878   kzE = (pzf[0]-vPosL[2])/(fXmiddle-vPosL[0]);
2879   //
2880   // get constrained parameters
2881   //
2882   Double_t xvertex=vPosL[0]-fXmiddle;
2883   fyf.AddPoint(&xvertex,vPosL[1], 0.00001);
2884   fzf.AddPoint(&xvertex,vPosL[2], 2.);
2885   fyf.Eval(); 
2886   fyf.GetParameters(pyfc); 
2887   fzf.Eval(); 
2888   fzf.GetParameters(pzfc); 
2889   //
2890   //
2891   // Make Fitters and params for 5 fitters
2892   // 1-4 OROC quadrants 
2893   //   0 IROC
2894   //
2895   static TLinearFitter *fittersY[5]={0,0,0,0,0};   // calgrind recomendation - fater to clear points
2896   static TLinearFitter *fittersZ[5]={0,0,0,0,0};   // than create the fitter
2897   if (fittersY[0]==0){
2898     for (Int_t i=0;i<5;i++) {
2899       fittersY[i] = new TLinearFitter(2,"pol1");
2900       fittersZ[i] = new TLinearFitter(2,"pol1");
2901     }
2902   }
2903   //
2904   Int_t npoints[5];
2905   TVectorD paramsY[5];
2906   TVectorD errorsY[5];
2907   TMatrixD covY[5];
2908   Double_t chi2CY[5];
2909   TVectorD paramsZ[5];
2910   TVectorD errorsZ[5];
2911   TMatrixD covZ[5];
2912   Double_t chi2CZ[5];
2913   for (Int_t i=0;i<5;i++) {
2914     npoints[i]=0;
2915     paramsY[i].ResizeTo(2);
2916     errorsY[i].ResizeTo(2);
2917     covY[i].ResizeTo(2,2);
2918     paramsZ[i].ResizeTo(2);
2919     errorsZ[i].ResizeTo(2);
2920     covZ[i].ResizeTo(2,2);
2921     fittersY[i]->ClearPoints();
2922     fittersZ[i]->ClearPoints();
2923   }
2924   //
2925   // Update fitters
2926   //
2927   Int_t countRes=0;
2928   for (Int_t irow=0;irow<159;irow++) {
2929     AliTPCclusterMI *c=track->GetClusterPointer(irow);
2930     if (!c) continue;
2931     if ((c->GetDetector()%36)!=isec) continue;
2932     Double_t dx = c->GetX()-fXmiddle;
2933     Double_t x[2]={dx, dx*dx};
2934     Double_t yfit  =  pyf[0]+pyf[1]*dx;
2935     Double_t zfit  =  pzf[0]+pzf[1]*dx;
2936     Double_t yfitC  =  pyfc[0]+pyfc[1]*dx;
2937     Double_t zfitC  =  pzfc[0]+pzfc[1]*dx;
2938     Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2939     if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
2940     if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
2941     if (dedge<kedgey) continue;
2942     Double_t corrtrY =  
2943       corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
2944                               c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2945     if (TMath::Abs(corrtrY)>kMaxCorrY) continue;  
2946     //
2947     vecX[countRes]=c->GetX();
2948     vecY[countRes]=c->GetY()-yfit;
2949     vecZ[countRes]=c->GetZ()-zfit;
2950     countRes++;
2951     //
2952     // Fill THnSparse cluster residuals
2953     // use only primary candidates with ITS signal
2954     if (fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){    
2955       Double_t resVector[5];
2956       resVector[1]= 9.*gphi/TMath::Pi();
2957       resVector[2]= TMath::Sqrt(c->GetX()*c->GetX()+c->GetY()*c->GetY());
2958       resVector[3]= c->GetZ()/resVector[2];
2959       //
2960       //
2961       resVector[0]= c->GetY()-yfitC;
2962       fClusterDelta[0]->Fill(resVector);
2963       resVector[0]= c->GetZ()-zfitC;
2964       fClusterDelta[1]->Fill(resVector);
2965     }
2966     if (c->GetRow()<kdrow0Fit) continue;      
2967     if (c->GetRow()>159-kdrow1Fit) continue;      
2968     //
2969
2970     if (c->GetDetector()>35){      
2971       if (c->GetX()<fXquadrant){
2972         if (yfit<-kPRFWidth) fittersY[1]->AddPoint(x,c->GetY(),0.1);
2973         if (yfit<-kPRFWidth) fittersZ[1]->AddPoint(x,c->GetZ(),0.1);
2974         if (yfit>kPRFWidth) fittersY[2]->AddPoint(x,c->GetY(),0.1);
2975         if (yfit>kPRFWidth) fittersZ[2]->AddPoint(x,c->GetZ(),0.1);
2976       }
2977       if (c->GetX()>fXquadrant){
2978         if (yfit<-kPRFWidth) fittersY[3]->AddPoint(x,c->GetY(),0.1);
2979         if (yfit<-kPRFWidth) fittersZ[3]->AddPoint(x,c->GetZ(),0.1);
2980         if (yfit>kPRFWidth) fittersY[4]->AddPoint(x,c->GetY(),0.1);
2981         if (yfit>kPRFWidth) fittersZ[4]->AddPoint(x,c->GetZ(),0.1);
2982       }
2983     }
2984     if (c->GetDetector()<36){
2985       fittersY[0]->AddPoint(x,c->GetY(),0.1);
2986       fittersZ[0]->AddPoint(x,c->GetZ(),0.1);
2987     }
2988   }
2989   //
2990   //
2991   //
2992   for (Int_t i=0;i<5;i++) {
2993     npoints[i] = fittersY[i]->GetNpoints();
2994     if (npoints[i]>=kMinClusterQ){
2995       // Y fit 
2996       fittersY[i]->Eval();
2997       Double_t chi2FacY = TMath::Sqrt(fittersY[i]->GetChisquare()/(fittersY[i]->GetNpoints()-2));
2998       chi2CY[i]=chi2FacY;
2999       fittersY[i]->GetParameters(paramsY[i]);
3000       fittersY[i]->GetErrors(errorsY[i]);
3001       fittersY[i]->GetCovarianceMatrix(covY[i]);
3002       // renormalize errors
3003       errorsY[i][0]*=chi2FacY;
3004       errorsY[i][1]*=chi2FacY;
3005       covY[i](0,0)*=chi2FacY*chi2FacY;
3006       covY[i](0,1)*=chi2FacY*chi2FacY;
3007       covY[i](1,0)*=chi2FacY*chi2FacY;
3008       covY[i](1,1)*=chi2FacY*chi2FacY;
3009       // Z fit
3010       fittersZ[i]->Eval();
3011       Double_t chi2FacZ = TMath::Sqrt(fittersZ[i]->GetChisquare()/(fittersZ[i]->GetNpoints()-2));
3012       chi2CZ[i]=chi2FacZ;
3013       fittersZ[i]->GetParameters(paramsZ[i]);
3014       fittersZ[i]->GetErrors(errorsZ[i]);
3015       fittersZ[i]->GetCovarianceMatrix(covZ[i]);
3016       // renormalize errors
3017       errorsZ[i][0]*=chi2FacZ;
3018       errorsZ[i][1]*=chi2FacZ;
3019       covZ[i](0,0)*=chi2FacZ*chi2FacZ;
3020       covZ[i](0,1)*=chi2FacZ*chi2FacZ;
3021       covZ[i](1,0)*=chi2FacZ*chi2FacZ;
3022       covZ[i](1,1)*=chi2FacZ*chi2FacZ;      
3023     }
3024   }
3025   //
3026   // void UpdateSectorKalman
3027   //
3028   for (Int_t i0=0;i0<5;i0++){
3029     for (Int_t i1=i0+1;i1<5;i1++){
3030       if(npoints[i0]<kMinClusterQ) continue;
3031       if(npoints[i1]<kMinClusterQ) continue;
3032       TMatrixD v0(4,1),v1(4,1);        // measurement
3033       TMatrixD cov0(4,4),cov1(4,4);    // covariance
3034       //
3035       v0(0,0)= paramsY[i0][0];         v1(0,0)= paramsY[i1][0]; 
3036       v0(1,0)= paramsY[i0][1];         v1(1,0)= paramsY[i1][1]; 
3037       v0(2,0)= paramsZ[i0][0];         v1(2,0)= paramsZ[i1][0]; 
3038       v0(3,0)= paramsZ[i0][1];         v1(3,0)= paramsZ[i1][1]; 
3039       //covariance i0
3040       cov0(0,0) = covY[i0](0,0);
3041       cov0(1,1) = covY[i0](1,1);
3042       cov0(1,0) = covY[i0](1,0);
3043       cov0(0,1) = covY[i0](0,1); 
3044       cov0(2,2) = covZ[i0](0,0);
3045       cov0(3,3) = covZ[i0](1,1);
3046       cov0(3,2) = covZ[i0](1,0);
3047       cov0(2,3) = covZ[i0](0,1);
3048       //covariance i1
3049       cov1(0,0) = covY[i1](0,0);
3050       cov1(1,1) = covY[i1](1,1);
3051       cov1(1,0) = covY[i1](1,0);
3052       cov1(0,1) = covY[i1](0,1); 
3053       cov1(2,2) = covZ[i1](0,0);
3054       cov1(3,3) = covZ[i1](1,1);
3055       cov1(3,2) = covZ[i1](1,0);
3056       cov1(2,3) = covZ[i1](0,1);
3057       //
3058       // And now update
3059       //
3060       if (TMath::Abs(pyf[1])<0.8){ //angular cut        
3061         UpdateSectorKalman(isec, i0,i1, &v0,&cov0,&v1,&cov1);
3062       }
3063     }
3064   }
3065
3066   //
3067   // Dump debug information
3068   //
3069   if (fStreamLevel>0){
3070     // get vertex position
3071      //
3072     TTreeSRedirector *cstream = GetDebugStreamer();  
3073
3074     
3075     if (cstream){      
3076       for (Int_t i0=0;i0<5;i0++){
3077         for (Int_t i1=i0;i1<5;i1++){
3078           if (i0==i1) continue;
3079           if(npoints[i0]<kMinClusterQ) continue;
3080           if(npoints[i1]<kMinClusterQ) continue;
3081           (*cstream)<<"sectorAlign"<<
3082             "run="<<fRun<<              //  run number
3083             "event="<<fEvent<<          //  event number
3084             "time="<<fTime<<            //  time stamp of event
3085             "trigger="<<fTrigger<<      //  trigger
3086             "triggerClass="<<&fTriggerClass<<      //  trigger
3087             "mag="<<fMagF<<             //  magnetic field        
3088             "isec="<<isec<<             //  current sector 
3089             //
3090             "xref="<<fXmiddle<<         // reference X
3091             "vPosG.="<<&vPosG<<        // vertex position in global system
3092             "vPosL.="<<&vPosL<<        // vertex position in local  system
3093             "vImpact.="<<&vImpact<<   // track impact parameter
3094             //"tofSignal="<<tofSignal<<   // tof signal
3095             "tpcPosG.="<<&tpcPosG<<   // global position of track at the middle of fXmiddle
3096             "lphi="<<lphi<<             // expected local phi angle - if vertex at 0
3097             "gphi="<<gphi<<             // expected global phi if vertex at 0
3098             "ky="<<ky<<
3099             "kyE="<<kyE<<               // expect ky - assiming pirmary track
3100             "kzE="<<kzE<<               // expected kz - assuming primary tracks
3101             "salpha="<<alpha<<          // sector alpha
3102             "scos="<<scos<<              // tracking cosinus
3103             "ssin="<<ssin<<              // tracking sinus
3104             //
3105             // full fit
3106             //
3107             "nf="<<nf<<                 //  total number of points
3108             "pyf.="<<&pyf<<             //  full OROC fit y
3109             "pzf.="<<&pzf<<             //  full OROC fit z
3110             "vX.="<<&vecX<<              //  x cluster
3111             "vY.="<<&vecY<<              //  y residual cluster
3112             "vZ.="<<&vecZ<<              //  z residual cluster
3113             // quadrant and IROC fit
3114             "i0="<<i0<<                 // quadrant number
3115             "i1="<<i1<<                 
3116             "n0="<<npoints[i0]<<        // number of points
3117             "n1="<<npoints[i1]<<
3118             //
3119             "py0.="<<&paramsY[i0]<<       // parameters
3120             "py1.="<<&paramsY[i1]<< 
3121             "ey0.="<<&errorsY[i0]<<       // errors
3122             "ey1.="<<&errorsY[i1]<<
3123             "chiy0="<<chi2CY[i0]<<       // chi2s
3124             "chiy1="<<chi2CY[i1]<<
3125             //
3126             "pz0.="<<&paramsZ[i0]<<       // parameters
3127             "pz1.="<<&paramsZ[i1]<< 
3128             "ez0.="<<&errorsZ[i0]<<       // errors
3129             "ez1.="<<&errorsZ[i1]<<
3130             "chiz0="<<chi2CZ[i0]<<       // chi2s
3131             "chiz1="<<chi2CZ[i1]<<
3132             "\n";
3133         }    
3134       }
3135     }
3136   }
3137 }
3138
3139 void AliTPCcalibAlign::UpdateSectorKalman(Int_t sector, Int_t quadrant0, Int_t quadrant1,  TMatrixD *const p0, TMatrixD *const c0, TMatrixD *const p1, TMatrixD *const c1 ){
3140   //
3141   //
3142   // tracks are refitted at sector middle
3143   //
3144   if (fArraySectorIntParam.At(0)==NULL) MakeSectorKalman();
3145   //
3146   //
3147   static TMatrixD matHk(4,30);    // vector to mesurement
3148   static TMatrixD measR(4,4);     // measurement error 
3149   //  static TMatrixD matQk(2,2);     // prediction noise vector
3150   static TMatrixD vecZk(4,1);     // measurement
3151   //
3152   static TMatrixD vecYk(4,1);     // Innovation or measurement residual
3153   static TMatrixD matHkT(30,4);   // helper matrix Hk transpose
3154   static TMatrixD matSk(4,4);     // Innovation (or residual) covariance
3155   static TMatrixD matKk(30,4);    // Optimal Kalman gain
3156   static TMatrixD mat1(30,30);    // update covariance matrix
3157   static TMatrixD covXk2(30,30);  // helper matrix
3158   //
3159   TMatrixD *vOrig = (TMatrixD*)(fArraySectorIntParam.At(sector));
3160   TMatrixD *cOrig = (TMatrixD*)(fArraySectorIntCovar.At(sector));
3161   //
3162   TMatrixD vecXk(*vOrig);             // X vector
3163   TMatrixD covXk(*cOrig);             // X covariance
3164   //
3165   //Unit matrix
3166   //
3167   for (Int_t i=0;i<30;i++)
3168     for (Int_t j=0;j<30;j++){
3169       mat1(i,j)=0;
3170       if (i==j) mat1(i,j)=1;
3171     }
3172   //
3173   //
3174   // matHk - vector to measurement
3175   //
3176   for (Int_t i=0;i<4;i++)
3177     for (Int_t j=0;j<30;j++){
3178       matHk(i,j)=0;
3179     }
3180   //
3181   // Measurement  
3182   // 0  - y
3183   // 1  - ky
3184   // 2  - z
3185   // 3  - kz
3186   
3187   matHk(0,6*quadrant1+4)  =  1.;            // delta y
3188   matHk(1,6*quadrant1+0)  =  1.;            // delta ky
3189   matHk(2,6*quadrant1+5)  =  1.;            // delta z
3190   matHk(3,6*quadrant1+1)  =  1.;            // delta kz
3191   // bug fix 24.02  - aware of sign in dx
3192   matHk(0,6*quadrant1+3)  =  -(*p0)(1,0);    // delta x to delta y  - through ky
3193   matHk(2,6*quadrant1+3)  =  -(*p0)(3,0);    // delta x to delta z  - thorugh kz
3194   matHk(2,6*quadrant1+2)  =  ((*p0)(0,0));  // y       to delta z  - through psiz
3195   //
3196   matHk(0,6*quadrant0+4)  =  -1.;           // delta y
3197   matHk(1,6*quadrant0+0)  =  -1.;           // delta ky
3198   matHk(2,6*quadrant0+5)  =  -1.;           // delta z
3199   matHk(3,6*quadrant0+1)  =  -1.;           // delta kz
3200   // bug fix 24.02 be aware of sign in dx
3201   matHk(0,6*quadrant0+3)  =  ((*p0)(1,0)); // delta x to delta y  - through ky
3202   matHk(2,6*quadrant0+3)  =  ((*p0)(3,0)); // delta x to delta z  - thorugh kz
3203   matHk(2,6*quadrant0+2)  =  -((*p0)(0,0)); // y       to delta z  - through psiz
3204
3205   //
3206   //
3207   
3208   vecZk =(*p1)-(*p0);                 // measurement
3209   measR =(*c1)+(*c0);                 // error of measurement
3210   vecYk = vecZk-matHk*vecXk;          // Innovation or measurement residual
3211   //
3212   //
3213   matHkT=matHk.T(); matHk.T();
3214   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
3215   matSk.Invert();
3216   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
3217   vecXk += matKk*vecYk;                    //  updated vector 
3218   covXk2= (mat1-(matKk*matHk));
3219   covXk =  covXk2*covXk;    
3220   //
3221   //
3222   (*cOrig)=covXk;
3223   (*vOrig)=vecXk;
3224 }
3225
3226 void AliTPCcalibAlign::MakeSectorKalman(){
3227   //
3228   // Make a initial Kalman paramaters for IROC - Quadrants alignment
3229   //
3230   TMatrixD param(5*6,1);
3231   TMatrixD covar(5*6,5*6);
3232   //
3233   // Set inital parameters
3234   //
3235   for (Int_t ip=0;ip<5*6;ip++) param(ip,0)=0;  // mean alignment to 0
3236   //
3237   for (Int_t iq=0;iq<5;iq++){
3238     // Initial uncertinty
3239     covar(iq*6+0,iq*6+0) = 0.002*0.002;        // 2 mrad  
3240     covar(iq*6+1,iq*6+1) = 0.002*0.002;        // 2 mrad  rotation
3241     covar(iq*6+2,iq*6+2) = 0.002*0.002;        // 2 mrad 
3242     //
3243     covar(iq*6+3,iq*6+3) = 0.02*0.02;        // 0.2 mm  
3244     covar(iq*6+4,iq*6+4) = 0.02*0.02;        // 0.2 mm  translation
3245     covar(iq*6+5,iq*6+5) = 0.02*0.02;        // 0.2 mm 
3246   }
3247
3248   for (Int_t isec=0;isec<36;isec++){
3249     fArraySectorIntParam.AddAt(param.Clone(),isec);
3250     fArraySectorIntCovar.AddAt(covar.Clone(),isec);
3251   }
3252 }
3253
3254 void AliTPCcalibAlign::UpdateSectorKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){
3255   //
3256   // Update kalman vector para0 with vector par1
3257   // Used for merging
3258   //
3259   static TMatrixD matHk(30,30);    // vector to mesurement
3260   static TMatrixD measR(30,30);     // measurement error 
3261   static TMatrixD vecZk(30,1);     // measurement
3262   //
3263   static TMatrixD vecYk(30,1);     // Innovation or measurement residual
3264   static TMatrixD matHkT(30,30);   // helper matrix Hk transpose
3265   static TMatrixD matSk(30,30);     // Innovation (or residual) covariance
3266   static TMatrixD matKk(30,30);    // Optimal Kalman gain
3267   static TMatrixD mat1(30,30);    // update covariance matrix
3268   static TMatrixD covXk2(30,30);  // helper matrix
3269   //
3270   TMatrixD vecXk(par0);             // X vector
3271   TMatrixD covXk(cov0);             // X covariance
3272
3273   //
3274   //Unit matrix
3275   //
3276   for (Int_t i=0;i<30;i++)
3277     for (Int_t j=0;j<30;j++){
3278       mat1(i,j)=0;
3279       if (i==j) mat1(i,j)=1;
3280     }
3281   matHk = mat1;                       // unit matrix 
3282   //
3283   vecZk = par1;                       // measurement
3284   measR = cov1;                        // error of measurement
3285   vecYk = vecZk-matHk*vecXk;          // Innovation or measurement residual
3286   //
3287   matHkT=matHk.T(); matHk.T();
3288   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
3289   matSk.Invert();
3290   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
3291   //matKk.Print();
3292   vecXk += matKk*vecYk;                    //  updated vector 
3293   covXk2= (mat1-(matKk*matHk));
3294   covXk =  covXk2*covXk;   
3295   CheckCovariance(covXk);
3296   CheckCovariance(cov1);
3297   //
3298   par0  = vecXk;                      // update measurement param
3299   cov0  = covXk;                      // update measurement covar
3300 }
3301
3302 Double_t AliTPCcalibAlign::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
3303   //
3304   // Get position correction for given sector
3305   //
3306
3307   TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
3308   if (!param) return 0;
3309   Int_t quadrant=0;
3310   if(lx>fXIO){
3311     if (lx<fXquadrant) {
3312       if (ly<0) quadrant=1;
3313       if (ly>0) quadrant=2;
3314     }
3315     if (lx>fXquadrant) {
3316       if (ly<0) quadrant=3;
3317       if (ly>0) quadrant=4;
3318     }
3319   }
3320   Double_t a10 = (*param)(quadrant*6+0,0);
3321   Double_t a20 = (*param)(quadrant*6+1,0);
3322   Double_t a21 = (*param)(quadrant*6+2,0);
3323   Double_t dx  = (*param)(quadrant*6+3,0);
3324   Double_t dy  = (*param)(quadrant*6+4,0);
3325   Double_t dz  = (*param)(quadrant*6+5,0);
3326   Double_t deltaX = lx-fXIO;
3327   if (coord==0) return dx;
3328   if (coord==1) return dy+deltaX*a10;
3329   if (coord==2) return dz+deltaX*a20+ly*a21;
3330   return 0;
3331 }
3332
3333 Double_t AliTPCcalibAlign::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
3334   //
3335   //
3336   //
3337   if (!Instance()) return 0;
3338   return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz);
3339 }
3340
3341 void AliTPCcalibAlign::MakeKalman(){
3342   //
3343   // Make a initial Kalman paramaters for sector Alignemnt
3344   //
3345   fSectorParamA = new TMatrixD(6*36+6,1);
3346   fSectorParamC = new TMatrixD(6*36+6,1);
3347   fSectorCovarA = new TMatrixD(6*36+6,6*36+6);
3348   fSectorCovarC = new TMatrixD(6*36+6,6*36+6);
3349   //
3350   // set starting parameters at 0
3351   //
3352   for (Int_t isec=0;isec<37;isec++)
3353     for (Int_t ipar=0;ipar<6;ipar++){
3354       (*fSectorParamA)(isec*6+ipar,0) =0;
3355       (*fSectorParamC)(isec*6+ipar,0) =0;
3356   }
3357   //
3358   // set starting covariance
3359   //
3360   for (Int_t isec=0;isec<36;isec++)
3361     for (Int_t ipar=0;ipar<6;ipar++){
3362       if (ipar<3){
3363         (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.002*0.002;   // 2 mrad
3364         (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.002*0.002;
3365       }
3366       if (ipar>=3){
3367         (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.02*0.02;     // 0.2 mm
3368         (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.02*0.02;
3369       }
3370   }
3371   (*fSectorCovarA)(36*6+0,36*6+0) =0.04;   // common shift y  up-up
3372   (*fSectorCovarA)(36*6+1,36*6+1) =0.04;   // common shift y  down-down
3373   (*fSectorCovarA)(36*6+2,36*6+2) =0.04;   // common shift y  up-down
3374   (*fSectorCovarA)(36*6+3,36*6+3) =0.004;   // common shift phi  up-up
3375   (*fSectorCovarA)(36*6+4,36*6+4) =0.004;   // common shift phi down-down
3376   (*fSectorCovarA)(36*6+5,36*6+5) =0.004;   // common shift phi up-down
3377   //
3378   (*fSectorCovarC)(36*6+0,36*6+0) =0.04;   // common shift y  up-up
3379   (*fSectorCovarC)(36*6+1,36*6+1) =0.04;   // common shift y  down-down
3380   (*fSectorCovarC)(36*6+2,36*6+2) =0.04;   // common shift y  up-down
3381   (*fSectorCovarC)(36*6+3,36*6+3) =0.004;   // common shift phi  up-up
3382   (*fSectorCovarC)(36*6+4,36*6+4) =0.004;   // common shift phi down-down
3383   (*fSectorCovarC)(36*6+5,36*6+5) =0.004;   // common shift phi up-down
3384 }
3385
3386 void AliTPCcalibAlign::UpdateKalman(Int_t sector0, Int_t sector1,  TMatrixD &p0, TMatrixD &c0, TMatrixD &p1, TMatrixD &c1){
3387   //
3388   // Update Kalman parameters
3389   // Note numbering from 0..36  0..17 IROC 18..35 OROC
3390   // 
3391   //
3392   if (fSectorParamA==NULL) MakeKalman();
3393   if (CheckCovariance(c0)>0) return;
3394   if (CheckCovariance(c1)>0) return;
3395   const Int_t nelem = 6*36+6;
3396   //
3397   //
3398   static TMatrixD matHk(4,nelem);    // vector to mesurement
3399   static TMatrixD measR(4,4);     // measurement error 
3400   static TMatrixD vecZk(4,1);     // measurement
3401   //
3402   static TMatrixD vecYk(4,1);     // Innovation or measurement residual
3403   static TMatrixD matHkT(nelem,4);   // helper matrix Hk transpose
3404   static TMatrixD matSk(4,4);     // Innovation (or residual) covariance
3405   static TMatrixD matKk(nelem,4);    // Optimal Kalman gain
3406   static TMatrixD mat1(nelem,nelem);    // update covariance matrix
3407   static TMatrixD covXk2(nelem,nelem);  // helper matrix
3408   //
3409   TMatrixD *vOrig = 0;
3410   TMatrixD *cOrig = 0;
3411   vOrig = (sector0%36>=18) ? fSectorParamA:fSectorParamC;
3412   cOrig = (sector0%36>=18) ? fSectorCovarA:fSectorCovarC;
3413   //
3414   Int_t sec0= sector0%18;
3415   Int_t sec1= sector1%18;
3416   if (sector0>35) sec0+=18;
3417   if (sector1>35) sec1+=18;
3418   //
3419   TMatrixD vecXk(*vOrig);             // X vector
3420   TMatrixD covXk(*cOrig);             // X covariance
3421   //
3422   //Unit matrix
3423   //
3424   for (Int_t i=0;i<nelem;i++)
3425     for (Int_t j=0;j<nelem;j++){
3426       mat1(i,j)=0;
3427       if (i==j) mat1(i,j)=1;
3428     }
3429   //
3430   //
3431   // matHk - vector to measurement
3432   //
3433   for (Int_t i=0;i<4;i++)
3434     for (Int_t j=0;j<nelem;j++){
3435       matHk(i,j)=0;
3436     }
3437   //
3438   // Measurement  
3439   // 0  - y
3440   // 1  - ky
3441   // 2  - z
3442   // 3  - kz
3443   
3444   matHk(0,6*sec1+4)  =  1.;            // delta y
3445   matHk(1,6*sec1+0)  =  1.;            // delta ky
3446   matHk(2,6*sec1+5)  =  1.;            // delta z
3447   matHk(3,6*sec1+1)  =  1.;            // delta kz
3448   matHk(0,6*sec1+3)  =  p0(1,0);    // delta x to delta y  - through ky
3449   matHk(2,6*sec1+3)  =  p0(3,0);    // delta x to delta z  - thorugh kz
3450   matHk(2,6*sec1+2)  =  p0(0,0);    // y       to delta z  - through psiz
3451   //
3452   matHk(0,6*sec0+4)  =  -1.;           // delta y
3453   matHk(1,6*sec0+0)  =  -1.;           // delta ky
3454   matHk(2,6*sec0+5)  =  -1.;           // delta z
3455   matHk(3,6*sec0+1)  =  -1.;           // delta kz
3456   matHk(0,6*sec0+3)  =  -p0(1,0); // delta x to delta y  - through ky
3457   matHk(2,6*sec0+3)  =  -p0(3,0); // delta x to delta z  - thorugh kz
3458   matHk(2,6*sec0+2)  =  -p0(0,0); // y       to delta z  - through psiz
3459
3460   Int_t dsec = (sector1%18)-(sector0%18);
3461   if (dsec<-2) dsec+=18; 
3462   if (TMath::Abs(dsec)==1){