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