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