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