Corrected UInt_t <-> Int_t conversion
[u/mrichter/AliRoot.git] / TPC / AliTPCkalmanAlign.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   TPC Kalman filter implementation for TPC sector alignment.
18   Output of the AliTPCcalibAlign is used as a input for TPC global alignment.
19   In AliTPCcalibAlign histograms - track parameter matching at sector boundaries are created.
20   Each sector is aligned with 5 neighborhoud (sectors)
21     1. Up-down    - 1
22     2. Left-right - 4
23
24   Sector alignment parameters are obtained finding the alignment parameters, minimizing
25   misalignmnet for all piars fo sectors.
26
27   Global minimization-  MakeGlobalAlign
28   
29
30   Example usage:
31   gSystem->Load("libANALYSIS");
32   gSystem->Load("libTPCcalib");
33   //
34   Int_t run=117220;
35   .x ConfigCalibTrain.C(run)
36   calibDB = AliTPCcalibDB::Instance()
37
38   AliTPCkalmanAlign kalmanAlign("TPC align", "TPC align");  // create the object
39   kalmanAlign.ReadAlign("CalibObjects.root");               // read the calibration form file         
40   kalmanAlign.MakeGlobalAlign();                            // do kalman alignment
41   kalmanAlign.DrawDeltaAlign();                             // make QA plot
42   //
43   
44
45 */
46 #include "TMath.h"
47 #include "TTreeStream.h"
48 #include "TGraph.h"
49 #include "TGraphErrors.h"
50 #include "TVectorD.h"
51 #include "TClonesArray.h"
52 #include "TCut.h"
53 #include "TChain.h"
54 #include "AliXRDPROOFtoolkit.h"
55 #include "TLegend.h"
56 #include "TCanvas.h"
57
58 #include "AliTPCcalibDB.h"
59 #include "AliTPCCalROC.h"
60 #include "AliCDBMetaData.h"
61 #include "AliCDBId.h"
62 #include "AliCDBManager.h"
63 #include "AliCDBStorage.h"
64 #include "AliCDBEntry.h"
65 #include "AliAlignObjParams.h"
66 #include "AliTPCROC.h"
67 #include "AliTracker.h"
68 #include "TFile.h"
69 #include "TLinearFitter.h"
70 #include "AliTPCcalibAlign.h"
71 #include "TH1.h"
72 #include "AliTPCCalPad.h"
73 #include "AliTPCkalmanAlign.h"
74 #include "TStatToolkit.h"
75 #include "AliTPCPreprocessorOnline.h"
76 #include "TPostScript.h"
77 #include "AliGRPObject.h"
78
79 AliTPCkalmanAlign::AliTPCkalmanAlign():
80   TNamed(),
81   fCalibAlign(0),     // kalman alignnmnt
82   fOriginalAlign(0),   // original alignment 0 read for the OCDB
83   fNewAlign(0),
84   fPadTime0(0),
85   fFitCEGlobal(0),
86   fFitCELocal(0)
87 {
88   //
89   // Default constructor
90   //
91   for (Int_t i=0; i<4; i++){
92     fDelta1D[i]=0;
93     fCovar1D[i]=0;
94   }
95 }
96
97 AliTPCkalmanAlign::AliTPCkalmanAlign(const char* name, const char* title): 
98   TNamed(name,title),
99   fCalibAlign(0),     // kalman alignnmnt
100   fOriginalAlign(0),   // original alignment 0 read for the OCDB
101   fNewAlign(0),
102   fPadTime0(0),
103   fFitCEGlobal(0),
104   fFitCELocal(0) 
105 {
106   //
107   // Default constructor
108   //
109   for (Int_t i=0; i<4; i++){
110     fDelta1D[i]=0;
111     fCovar1D[i]=0;
112   }
113   fFitCEGlobal = new TObjArray(6); 
114   fFitCELocal  = new TObjArray(6); 
115   for (Int_t ipar=0; ipar<6;ipar++){
116     fFitCEGlobal->AddAt(new TVectorD(36),ipar);
117     fFitCELocal->AddAt(new TVectorD(36),ipar);
118   }
119 }
120
121 void AliTPCkalmanAlign::ReadAlign(const char *fname){
122   //
123   // Read old alignment used in the reco
124   // and the residual histograms
125   // WE ASSUME that the OCDB path is set in the same way as done in the calibration
126   //
127   TFile fcalib(fname);
128   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
129   fCalibAlign=0;
130   if (array) fCalibAlign=( AliTPCcalibAlign *)array->FindObject("alignTPC");
131   fCalibAlign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
132   //
133   // old alignment used
134   AliCDBEntry * cdbEntry= AliCDBManager::Instance()->Get("TPC/Align/Data");
135   fOriginalAlign =0;
136   if (cdbEntry){
137     fOriginalAlign = (TClonesArray*)cdbEntry->GetObject();
138   }
139   
140 }
141
142 void AliTPCkalmanAlign::BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t mean, Double_t sigma){
143   //
144   // Book Align 1D parameters and covariance
145   //
146   param.ResizeTo(72,1);
147   covar.ResizeTo(72,72);
148   for (Int_t i=0;i<72;i++){
149     param(i,0)=mean;
150     for (Int_t j=0;j<72;j++) covar(i,j)=0;      
151     covar(i,i)=sigma*sigma;
152   }
153 }
154
155
156 void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &vecXk, TMatrixD &covXk){
157   //
158   // Update 1D kalman filter
159   //
160   const Int_t knMeas=1;
161   const Int_t knElem=72;
162   TMatrixD mat1(72,72);            // update covariance matrix
163   TMatrixD matHk(1,knElem);        // vector to mesurement
164   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
165   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
166   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
167   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
168   TMatrixD covXk2(knElem,knElem);  // helper matrix
169   TMatrixD covXk3(knElem,knElem);  // helper matrix
170   TMatrixD vecZk(1,1);
171   TMatrixD measR(1,1);
172   vecZk(0,0)=delta;
173   measR(0,0)=sigma*sigma;
174   //
175   // reset matHk
176   for (Int_t iel=0;iel<knElem;iel++) 
177     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
178   //mat1
179   for (Int_t iel=0;iel<knElem;iel++) {
180     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
181     mat1(iel,iel)=1;
182   }
183   //
184   matHk(0, s1)=1;
185   matHk(0, s2)=-1;
186   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
187   matHkT=matHk.T(); matHk.T();
188   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
189   matSk.Invert();
190   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
191   vecXk += matKk*vecYk;                    //  updated vector 
192   covXk2= (mat1-(matKk*matHk));
193   covXk3 =  covXk2*covXk;          
194   covXk = covXk3;  
195   Int_t nrows=covXk3.GetNrows();
196   
197   for (Int_t irow=0; irow<nrows; irow++)
198     for (Int_t icol=0; icol<nrows; icol++){
199       // rounding problems - make matrix again symteric
200       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; 
201     }
202 }
203
204 void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
205   //
206   // Update 1D kalman filter - with one measurement
207   //
208   const Int_t knMeas=1;
209   const Int_t knElem=72;
210   TMatrixD mat1(72,72);            // update covariance matrix
211   TMatrixD matHk(1,knElem);        // vector to mesurement
212   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
213   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
214   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
215   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
216   TMatrixD covXk2(knElem,knElem);  // helper matrix
217   TMatrixD covXk3(knElem,knElem);  // helper matrix
218   TMatrixD vecZk(1,1);
219   TMatrixD measR(1,1);
220   vecZk(0,0)=delta;
221   measR(0,0)=sigma*sigma;
222   //
223   // reset matHk
224   for (Int_t iel=0;iel<knElem;iel++) 
225     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
226   //mat1
227   for (Int_t iel=0;iel<knElem;iel++) {
228     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
229     mat1(iel,iel)=1;
230   }
231   //
232   matHk(0, s1)=1;
233   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
234   matHkT=matHk.T(); matHk.T();
235   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
236   matSk.Invert();
237   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
238   vecXk += matKk*vecYk;                    //  updated vector 
239   covXk2= (mat1-(matKk*matHk));
240   covXk3 =  covXk2*covXk;          
241   covXk = covXk3;  
242   Int_t nrows=covXk3.GetNrows();
243   
244   for (Int_t irow=0; irow<nrows; irow++)
245     for (Int_t icol=0; icol<nrows; icol++){
246       // rounding problems - make matrix again symteric
247       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; 
248     }
249 }
250
251
252
253
254 void AliTPCkalmanAlign::MakeGlobalAlign(){
255   //
256   // Combine all pairs of fitters and make global alignemnt
257   //
258   
259   AliTPCkalmanAlign &kalmanAlign=*this;
260   TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
261   Int_t run = AliCDBManager::Instance()->GetRun();
262   AliGRPObject * grp = AliTPCcalibDB::Instance()->GetGRP(run);
263   Float_t bz = AliTracker::GetBz();
264   UInt_t timeS = grp->GetTimeStart();
265   UInt_t timeE = grp->GetTimeEnd();
266   UInt_t  time = (timeS+timeE)/2;
267   
268   //
269   // get ce info
270   //
271   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
272   TVectorD paramCE[72];
273   TMatrixD covarCE[72];
274   Int_t  statUpDown=0;   // statistic up down
275   Int_t  statLeftRight=0;   // statistic left-right
276   Float_t chi2;
277   for (Int_t isec=0; isec<72; isec++){
278     AliTPCCalROC * roc0  = padTime0->GetCalROC(isec);
279     roc0->GlobalFit(0,kFALSE,paramCE[isec],covarCE[isec],chi2,0);
280     (*pcstream)<<"ceFit"<<
281       "isec="<<isec<<
282       "p0.="<<&paramCE[isec]<<
283       "\n";
284   }
285
286   DumpOldAlignment(pcstream);
287   const Int_t kMinEntries=400;
288   TMatrixD vec[5];
289   TMatrixD cov[5];
290   kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
291   kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
292   kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
293   kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
294   //
295   TVectorD delta(5);
296   TVectorD rms(5);
297   TVectorD stat(5);
298   TH1 * his=0;
299   for (Int_t is0=0;is0<72;is0++)
300     for (Int_t is1=0;is1<72;is1++){
301       //
302       //
303       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
304       if (!his) continue;
305       if (his->GetEntries()<kMinEntries) continue;
306       delta[0]=his->GetMean();
307       rms[0]=his->GetRMS();
308       stat[0]=his->GetEntries();
309       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
310       //     
311       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
312       if (!his) continue;
313       delta[1]=his->GetMean();
314       rms[1]=his->GetRMS();
315       stat[1]=his->GetEntries();
316       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
317       //
318       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
319       if (!his) continue;
320       delta[2] = his->GetMean();
321       rms[2]=his->GetRMS();
322       stat[2]=his->GetEntries();
323       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
324       //
325       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
326       if (!his) continue;
327       delta[3] = his->GetMean();       
328       rms[3]=his->GetRMS();
329       stat[3]=his->GetEntries();
330       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
331       if (is1==is0+36) statUpDown+=Int_t(stat[0]);
332       if (is1==is0+35) statLeftRight+=Int_t(stat[0]);
333     }  
334
335   fDelta1D[0] = (TMatrixD*)vec[0].Clone();
336   fDelta1D[1] = (TMatrixD*)vec[1].Clone();
337   fDelta1D[2] = (TMatrixD*)vec[2].Clone();
338   fDelta1D[3] = (TMatrixD*)vec[3].Clone();
339   //
340   fCovar1D[0] = (TMatrixD*)cov[0].Clone();
341   fCovar1D[1] = (TMatrixD*)cov[1].Clone();
342   fCovar1D[2] = (TMatrixD*)cov[2].Clone();
343   fCovar1D[3] = (TMatrixD*)cov[3].Clone();
344   statUpDown/=36;
345   statLeftRight/=36;
346   MakeNewAlignment(kTRUE);
347   FitCE();
348   for (Int_t is0=0;is0<72;is0++)
349     for (Int_t is1=0;is1<72;is1++){
350       Bool_t isPair=kFALSE;
351       if (TMath::Abs(is0%18-is1%18)<2) isPair=kTRUE;
352       if (TMath::Abs(is0%18-is1%18)==17) isPair=kTRUE;
353       if (!isPair) continue;
354       stat[0]=0; stat[1]=0; stat[2]=0; stat[3]=0; 
355       //
356       //
357       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
358       if (his){
359         delta[0]=his->GetMean();
360         rms[0]=his->GetRMS();
361         stat[0]=his->GetEntries();
362       }
363       //     
364       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
365       if (his) {
366         delta[1]=his->GetMean();
367         rms[1]=his->GetRMS();
368         stat[1]=his->GetEntries();
369       }
370       //
371       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
372       if (his){
373         delta[2] = his->GetMean();
374         rms[2]=his->GetRMS();
375         stat[2]=his->GetEntries();
376       }
377       //
378       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
379       if (his){
380         delta[3] = his->GetMean();       
381         rms[3]=his->GetRMS();
382         stat[3]=his->GetEntries();
383       }
384       TVectorD fceG[8],fceL[6];
385       if (fFitCEGlobal)
386         for (Int_t ipar=0; ipar<8;ipar++){
387           fceG[ipar].ResizeTo(36);
388           if (ipar<6) fceL[ipar].ResizeTo(36);
389           if (fFitCEGlobal->At(ipar)){
390             fceG[ipar]=*((TVectorD*)fFitCEGlobal->At(ipar));
391             if (ipar<6){
392               fceL[ipar]=*((TVectorD*)fFitCELocal->At(ipar));
393             }
394           }           
395         }
396       
397       (*pcstream)<<"kalmanAlignDebug"<<
398         "run="<<run<<            // runs
399         "time="<<time<<          // time
400         "timeE="<<timeE<<        // sart of tun time
401         "timeS="<<timeS<<        // end od run time
402         "bz="<<bz<<
403         "is0="<<is0<<
404         "is1="<<is1<<
405         "delta.="<<&delta<<      // alignment deltas
406         "rms.="<<&rms<<          // rms
407         "stat.="<<&stat<<
408         "vec0.="<<&vec[0]<<
409         "vec1.="<<&vec[1]<<
410         "vec2.="<<&vec[2]<<
411         "vec3.="<<&vec[3]<<
412         "pceIn0.="<<&paramCE[is0%36]<<        // default CE parameters
413         "pceOut0.="<<&paramCE[is0%36+36]<<
414         "pceIn1.="<<&paramCE[is1%36]<<
415         "pceOut1.="<<&paramCE[is1%36+36]<<
416         //                     // current CE parameters form last calibration - not used in Reco
417         "fceG0.="<<&fceG[0]<<  // global fit of CE  - offset
418         "fceG1.="<<&fceG[1]<<  // global fit of CE  - Gy gradient 
419         "fceG2.="<<&fceG[2]<<  // global fit of CE  - Gx gradient
420         "fceG3.="<<&fceG[3]<<  // global fit of CE  - IROC-OROC offset
421         "fceG4.="<<&fceG[4]<<  // global fit of CE  - commont slope LX
422         "fceG5.="<<&fceG[5]<<  // global fit of CE  - delta slope LX
423         "fceG6.="<<&fceG[6]<<  // global fit of CE  - common slope LY
424         "fceG7.="<<&fceG[7]<<  // global fit of CE  - delta slope LY
425         //
426         "fceL0.="<<&fceL[0]<<  // local fit of CE  - offset to mean
427         "fceL1.="<<&fceL[1]<<  // local fit of CE  - IROC-OROC offset
428         "fceL2.="<<&fceL[2]<<  // local fit of CE  - common slope LX
429         "fceL3.="<<&fceL[3]<<  // local fit of CE  - delta slope  LX
430         "fceL4.="<<&fceL[4]<<  // local fit of CE  - common slope LY
431         "fceL5.="<<&fceL[5]<<  // local fit of CE  - delta slope  LY
432         "\n";
433     }
434   
435   (*pcstream)<<"runSummary"<<
436     "run="<<run<<                      // run number 
437     "bz="<<bz<<                        // bz field
438     "statUpDown="<<statUpDown<<        // stat up-down
439     "statLeftRight="<<statLeftRight<<  // stat left-right
440     "\n";
441
442   delete pcstream;
443 }
444
445
446
447
448
449
450 void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad  *pad, Int_t ustartRun, Int_t uendRun,  const char* storagePath ){
451   //
452   // Update OCDB
453   // .x ConfigCalibTrain.C(117116)
454   // AliTPCcalibDB::Instance()->GetPulserTmean()
455   // pad->Add(-pad->GetMean())
456   AliCDBMetaData *metaData= new AliCDBMetaData();
457   metaData->SetObjectClassName("TObjArray");
458   metaData->SetResponsible("Marian Ivanov");
459   metaData->SetBeamPeriod(1);
460   metaData->SetAliRootVersion("05-25-01"); //root version
461   metaData->SetComment("Calibration of the z - time misalignment");
462   AliCDBId* id1=NULL;
463   id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
464   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
465   gStorage->Put(pad, (*id1), metaData);
466 }
467
468
469
470 void AliTPCkalmanAlign::DrawDeltaAlign(){
471   //
472   // Draw the reuslts of the alignment
473   // Residual misalignment in respect with previous alignment shown
474   //
475   //
476   TFile f("AliTPCkalmanAlign.root","update");
477   TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
478   TH1::AddDirectory(0);
479   // 
480   treeDelta->SetAlias("sector","is0");
481   treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
482   treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
483   treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
484   treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
485   //
486   treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
487   treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
488   treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
489   treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
490   //
491   treeDelta->SetMarkerStyle(25);
492   treeDelta->SetMarkerColor(4);
493   treeDelta->SetLineColor(4);
494   const char *type[3]={"up-down","left-right","right-left"};
495   const char *gropt[3]={"alp","lp","lp"};
496   const char *varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10*vec0.fElements[is0]:sector"};
497   const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
498   const char *varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10*vec1.fElements[is0]:sector"};
499   const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
500   const char *varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000*vec3.fElements[is0]:sector"};
501   const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
502   const char *varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000*vec2.fElements[is0]:sector"};
503   const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
504   TLegend *legend = 0;
505   Int_t grcol[3]={2,1,4};
506   Int_t entries=0;
507   TGraph *grDelta[3]={0,0,0};
508   TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
509   TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
510   TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
511   TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
512   canvasDy->Divide(2,2);
513   canvasDz->Divide(2,2);
514   canvasDtheta->Divide(2,2);
515   canvasDphi->Divide(2,2);
516
517   //
518   // Dy
519   //
520   canvasDy->cd(1);
521   treeDelta->Draw("dYmeas:dYfit");
522   for (Int_t itype=0; itype<3; itype++){
523     canvasDy->cd(itype+2);
524     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
525     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
526     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
527     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
528     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
529     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
530     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
531     for (Int_t i=0; i<3; i++) {
532       grDelta[i]->SetMaximum(1.5);
533       grDelta[i]->SetMinimum(-1);
534       grDelta[i]->SetTitle(type[i]);
535       grDelta[i]->SetMarkerColor(grcol[i]); 
536       grDelta[i]->SetLineColor(grcol[i]); 
537       grDelta[i]->SetMarkerStyle(25+i); 
538       grDelta[i]->GetXaxis()->SetTitle("sector"); 
539       grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)"); 
540       if (itype==2 && i>0) continue;
541       grDelta[i]->Draw(gropt[i]); 
542       legend->AddEntry(grDelta[i]);
543     }
544     legend->Draw();
545   }
546   //
547   // Dz
548   //
549   canvasDz->cd();
550   canvasDz->cd(1);
551   treeDelta->Draw("dZmeas:dZfit");
552   for (Int_t itype=0; itype<3; itype++){
553     canvasDz->cd(itype+2);
554     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
555     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
556     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
557     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
558     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
559     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
560     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
561     for (Int_t i=0; i<3; i++) {
562       grDelta[i]->SetMaximum(1.5);
563       grDelta[i]->SetMinimum(-1);
564       grDelta[i]->SetTitle(type[i]);
565       grDelta[i]->SetMarkerColor(grcol[i]); 
566       grDelta[i]->SetLineColor(grcol[i]); 
567       grDelta[i]->SetMarkerStyle(25+i); 
568       grDelta[i]->GetXaxis()->SetTitle("sector"); 
569       grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)"); 
570       if (itype==2 && i>0) continue;
571       grDelta[i]->Draw(gropt[i]); 
572       legend->AddEntry(grDelta[i]);
573     }
574     legend->Draw();
575   }
576
577   //
578   // Dtheta
579   //
580   canvasDtheta->cd(1);
581   treeDelta->Draw("dThetameas:dThetafit");
582   for (Int_t itype=0; itype<3; itype++){
583     canvasDtheta->cd(itype+2);
584     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
585     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
586     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
587     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
588     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
589     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
590     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
591     for (Int_t i=0; i<3; i++) {
592       grDelta[i]->SetMaximum(4.);
593       grDelta[i]->SetMinimum(-3.);
594       grDelta[i]->SetTitle(type[i]);
595       grDelta[i]->SetMarkerColor(grcol[i]); 
596       grDelta[i]->SetLineColor(grcol[i]); 
597       grDelta[i]->SetMarkerStyle(25+i); 
598       grDelta[i]->GetXaxis()->SetTitle("sector"); 
599       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)"); 
600       if (itype==2 && i>0) continue;
601       grDelta[i]->Draw(gropt[i]); 
602       legend->AddEntry(grDelta[i]);
603     }
604     legend->Draw();
605   }
606
607   //
608   // Dphi
609   //
610   canvasDphi->cd();
611   canvasDphi->cd(1);
612   treeDelta->Draw("dPhimeas:dPhifit");
613   for (Int_t itype=0; itype<3; itype++){
614     canvasDphi->cd(itype+2);
615     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
616     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
617     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
618     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
619     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
620     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
621     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
622     for (Int_t i=0; i<3; i++) {
623       grDelta[i]->SetMaximum(4.);
624       grDelta[i]->SetMinimum(-3.);
625       grDelta[i]->SetTitle(type[i]);
626       grDelta[i]->SetMarkerColor(grcol[i]); 
627       grDelta[i]->SetLineColor(grcol[i]); 
628       grDelta[i]->SetMarkerStyle(25+i); 
629       grDelta[i]->GetXaxis()->SetTitle("sector"); 
630       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)"); 
631       if (itype==2 && i>0) continue;
632       grDelta[i]->Draw(gropt[i]); 
633       legend->AddEntry(grDelta[i]);
634     }
635     legend->Draw();
636   }
637   //
638   //
639   f.cd();
640   canvasDy->Write();
641   canvasDz->Write();
642   canvasDtheta->Write();
643   canvasDphi->Write();
644   //  
645   //
646   //
647   TPostScript *ps = new TPostScript("alignReport.ps", 112); 
648   ps->NewPage();
649   canvasDy->Draw();
650   ps->NewPage();
651   canvasDz->Draw();
652   ps->NewPage();
653   canvasDtheta->Draw();
654   ps->NewPage();
655   canvasDphi->Draw();
656   ps->Close();
657   delete ps;
658 }
659
660
661
662 void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
663   // Dump the content of old alignemnt
664   // Expected that the old alignmnet is loaded
665   //
666   if (!fOriginalAlign) return;
667   //
668   TVectorD localTrans(3);
669   TVectorD globalTrans(3);
670   TVectorD localRot(3);
671   TVectorD globalRot(3);
672   AliGeomManager::ELayerID idLayer;
673   Int_t idModule=0;
674   //
675   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
676     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
677     params->GetVolUID(idLayer,idModule);
678     params->GetLocalTranslation(localTrans.GetMatrixArray());
679     params->GetLocalAngles(localRot.GetMatrixArray());
680     params->GetTranslation(globalTrans.GetMatrixArray());
681     params->GetAngles(globalRot.GetMatrixArray());
682     Int_t sector=idModule;
683     if (idLayer>7) sector+=36;
684     (*pcstream)<<"oldAlign"<<
685       //"idLayer="<<idLayer<<
686       "idModule="<<idModule<<
687       "sector="<<sector<<
688       "lT.="<<&localTrans<<
689       "gT.="<<&localTrans<<
690       "lR.="<<&localRot<<
691       "gR.="<<&globalRot<<
692       "\n";
693   }
694 }
695
696
697 void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
698   //
699   // make a new Alignment entry
700   //
701   if (!fOriginalAlign) return;
702   //
703   TVectorD localTrans(3);
704   TVectorD globalTrans(3);
705   TVectorD localRot(3);
706   TVectorD globalRot(3);
707   //
708   TVectorD localTransNew(3);   // new entries
709   TVectorD globalTransNew(3);
710   TVectorD localRotNew(3);
711   TVectorD globalRotNew(3);
712   //
713   AliGeomManager::ELayerID idLayer;
714   Int_t idModule=0;
715   //
716   fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
717   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
718     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
719     //AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
720     params->GetVolUID(idLayer,idModule);
721     Int_t sector=(Int_t)idModule;
722     if (idLayer>7) sector+=36;
723     params->GetLocalTranslation(localTrans.GetMatrixArray());
724     params->GetLocalAngles(localRot.GetMatrixArray());
725     params->GetTranslation(globalTrans.GetMatrixArray());
726     params->GetAngles(globalRot.GetMatrixArray());
727     //
728     //
729     //
730     if (badd){ // addition if 
731       localTransNew=localTrans;
732       localRotNew=localRot;
733     }
734     localTransNew[1]=localTransNew[1]-((*fDelta1D[0])(sector,0));
735     localRot[0]  =localRot[0]-(*fDelta1D[2])(sector,0);
736     //
737     if (pcstream) (*pcstream)<<"alignParams"<<
738       //"idLayer="<<idLayer<<
739       "idModule="<<idModule<<
740       "sector="<<sector<<
741       "olT.="<<&localTrans<<
742       "olR.="<<&localRot<<
743       "ogT.="<<&localTrans<<
744       "ogR.="<<&globalRot<<
745       "nlT.="<<&localTransNew<<
746       "nlR.="<<&localRotNew<<
747       "ngT.="<<&localTransNew<<
748       "ngR.="<<&globalRotNew<<
749       "\n";
750   }
751 }
752
753
754
755 void AliTPCkalmanAlign::DrawAlignmentTrends(){
756   //
757   // Draw trends of alingment variables
758   //
759   /*
760     //1.  Create a list of  align data
761     //
762     //2. Filter list
763     AliXRDPROOFtoolkit::FilterList("align.list","AliTPCkalmanAlign.root kalmanAlignDebug",0);
764
765   */
766   AliXRDPROOFtoolkit toolkit;
767   TChain * chain = toolkit.MakeChainRandom("align.list.Good","kalmanAlignDebug",0,2000);
768   TChain * chainRef = toolkit.MakeChainRandom("alignRef.list","kalmanAlignDebug",0,2000);
769   chain->AddFriend(chainRef,"R");
770   chainRef->AddFriend(chainRef,"T");
771   //cuts
772   TCut cutS="stat.fElements[0]>200&&stat.fElements[1]>200&&stat.fElements[3]>200&&stat.fElements[3]>200";   //statistic in the bin
773   TCut cutST="T.stat.fElements[0]>200&&T.stat.fElements[1]>200&&T.stat.fElements[3]>200&&T.stat.fElements[3]>200";   //statistic in the bin
774   //  TTree *tree  = chain->CopyTree(cutS);
775   //TTree *treeR = chainRef->CopyTree(cutST);
776
777   TCanvas * canvasDy= new TCanvas("canvasDy","canvasDy");
778   TH1 *his=0;
779   TLegend *legend = 0;
780   //  Int_t grcol[3]={2,1,4};
781   
782   legend = new TLegend(0.7,0.6,0.9,0.9, "Alignment #Delta_{y}- Up-Down");
783   for (Int_t isec=0; isec<18; isec+=2){
784     chain->SetMarkerColor(1+(isec%5));
785     chain->SetMarkerStyle(isec+20);
786     chain->Draw("10*(delta.fElements[0]-R.delta.fElements[0]):run",cutS+Form("is1==is0+36&&is0==%d",isec),"profgoff");
787     his = (TH1*)(chain->GetHistogram()->Clone());
788     his->SetName(Form("#Delta_{Y} sector %d",isec));
789     his->SetTitle(Form("#Delta_{Y} sector %d",isec));
790     his->SetMaximum(1.);
791     his->SetMinimum(-1.);
792     his->GetYaxis()->SetTitle("#Delta_{y} (mm)");
793     his->GetXaxis()->SetTitle("run Number");
794     if (isec==0) his->Draw("");
795     if (isec>0) his->Draw("same");
796     legend->AddEntry(his);
797   }
798   legend->Draw();
799   canvasDy->Draw();
800 }
801
802
803
804
805
806
807 void AliTPCkalmanAlign::FitCE(){
808   //
809   // fit CE 
810   // 1. Global fit - gy and gx
811   // 2. Local X fit common 
812   // 3. Sector fit 
813   //
814   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
815   //
816   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
817   AliTPCCalPad *padNoise = AliTPCcalibDB::Instance()->GetPadNoise();
818   AliTPCCalPad * ceTmean = AliTPCcalibDB::Instance()->GetCETmean();  // CE information
819   AliTPCCalPad * ceTrms  = AliTPCcalibDB::Instance()->GetCETrms();
820   AliTPCCalPad * ceQmean = AliTPCcalibDB::Instance()->GetCEQmean();  
821   AliTPCCalPad * pulserTmean = AliTPCcalibDB::Instance()->GetPulserTmean(); //
822   AliTPCCalPad * pulserTrms = AliTPCcalibDB::Instance()->GetPulserTrms();
823   AliTPCCalPad * pulserQmean = AliTPCcalibDB::Instance()->GetPulserQmean();
824   AliTPCCalPad * dmap0  = AliTPCcalibDB::Instance()->GetDistortionMap(0);   // distortion maps
825   AliTPCCalPad * dmap1  = AliTPCcalibDB::Instance()->GetDistortionMap(1);
826   AliTPCCalPad * dmap2  = AliTPCcalibDB::Instance()->GetDistortionMap(2);
827   pulserTmean->Add(-pulserTmean->GetMean());
828   //
829   preprocesor->AddComponent(padTime0->Clone());
830   preprocesor->AddComponent(padNoise->Clone());
831   preprocesor->AddComponent(pulserTmean->Clone());
832   preprocesor->AddComponent(pulserQmean->Clone());
833   preprocesor->AddComponent(pulserTrms->Clone());
834   preprocesor->AddComponent(ceTmean->Clone());
835   preprocesor->AddComponent(ceQmean->Clone());
836   preprocesor->AddComponent(ceTrms->Clone());
837   preprocesor->AddComponent(dmap0->Clone());
838   preprocesor->AddComponent(dmap1->Clone());
839   preprocesor->AddComponent(dmap2->Clone());
840   preprocesor->DumpToFile("cetmean.root");
841
842   TCut cutNoise="abs(PadNoise.fElements/PadNoise_Median-1)<0.3";
843   TCut cutPulserT="abs(PulserTrms.fElements/PulserTrms_Median-1)<0.2";
844   TCut cutPulserQ="abs(PulserQmean.fElements/PulserQmean_Median-1)<0.2";
845   TCut cutCEQ="CEQmean.fElements>50";
846   TCut cutCET="abs(CETmean.fElements)<2";
847   TCut cutAll=cutNoise+cutPulserT+cutPulserQ+cutCEQ+cutCET;
848   //
849   //
850   TFile * f = new TFile("cetmean.root");
851   TTree * chain = (TTree*) f->Get("calPads");
852   Int_t entries = chain->Draw("1",cutAll,"goff");
853   if (entries<200000) return;  // no calibration available - pulser or CE or noise
854
855   TStatToolkit toolkit;
856   Double_t chi2=0;
857   Int_t    npoints=0;
858   TVectorD param;
859   TMatrixD covar;
860   //
861   // make a aliases
862   AliTPCkalmanAlign::MakeAliasCE(chain);
863   TString  fstringG="";              // global part
864   //
865   fstringG+="Gy++";                  // 1 - global y
866   fstringG+="Gx++";                  // 2 - global x
867   // 
868   fstringG+="isin++";                // 3 -delta IROC-OROC offset
869   fstringG+="Lx++";                  // 4 -common slope 
870   fstringG+="Lx*isin++";             // 5 -delta slope 
871   fstringG+="Ly++";                  // 6 -common slope 
872   fstringG+="Ly*isin++";             // 7 -delta slope 
873   TVectorD vecG[2];
874   TString * strFitG=0;
875   TString * strFitLX=0;
876   //
877   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideA"+cutAll, chi2,npoints,vecG[0],covar,-1,0, 10000000, kFALSE);
878   chain->SetAlias("tfitGA",strFitG->Data());
879   strFitG->Tokenize("++")->Print();
880   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
881   //
882   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideC"+cutAll, chi2,npoints,vecG[1],covar,-1,0, 10000000, kFALSE);
883   chain->SetAlias("tfitGC",strFitG->Data());
884   strFitG->Tokenize("++")->Print();
885   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
886   //
887   AliTPCCalPad *padFitG =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[0],vecG[1]);
888   AliTPCCalPad *padFitLX=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[0],vecG[1]);
889   // swap a side and c side
890   AliTPCCalPad *padFitGSwap =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[1],vecG[0]);
891   AliTPCCalPad *padFitLXSwap=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[1],vecG[0]);
892   padFitG->SetName("CEG");
893   padFitLX->SetName("CELX");
894   padFitGSwap->SetName("CEGS");
895   padFitLXSwap->SetName("CELXS");
896   preprocesor->AddComponent(padFitG->Clone());
897   preprocesor->AddComponent(padFitLX->Clone());
898   preprocesor->AddComponent(padFitGSwap->Clone());
899   preprocesor->AddComponent(padFitLXSwap->Clone());
900   preprocesor->DumpToFile("cetmean.root");   // add it to the file
901   //
902   // make local fits
903   //
904   f = new TFile("cetmean.root");
905   chain = (TTree*) f->Get("calPads");
906   AliTPCkalmanAlign::MakeAliasCE(chain);
907   TString  fstringL="";              // local fit
908   //                                 // 0. delta common
909   fstringL+="isin++";                // 1. delta IROC-OROC offset
910   fstringL+="Lx++";                  // 2. common slope 
911   fstringL+="Lx*isin++";             // 3. delta slope 
912   fstringL+="Ly++";                  // 4. common slope 
913   fstringL+="Ly*isin++";             // 5. delta slope 
914   TVectorD vecL[36];
915   TVectorD dummy(6);
916   AliTPCCalPad *padFitLCE = new AliTPCCalPad("LocalCE","LocalCE");
917   AliTPCCalPad *padFitTmpCE;
918   for (Int_t isec=0; isec<36; isec++){
919     TCut cutSector=Form("(sector%%36)==%d",isec);
920     strFitLX = TStatToolkit::FitPlane(chain,"deltaT-CEG.fElements-CELX.fElements", fstringL.Data(),cutSector+cutAll+"abs(deltaT-CEG.fElements-CELX.fElements)<0.4", chi2,npoints,vecL[isec],covar,-1,0, 10000000, kFALSE);
921     printf("sec=%d\tchi2=%f\n",isec,TMath::Sqrt(chi2/npoints));
922     delete strFitLX;
923     //
924     TString fitL=Form("((sector%%36)==%d)++((sector%%36)==%d)*(sector<36)++((sector%%36)==%d)*(lx-133)/100.++((sector%%36)==%d)*(sector<36)*(lx-133)/100.++((sector%%36)==%d)*(ly)/100.++((sector%%36)==%d)*(sector<36)*(ly)/100.",isec,isec,isec,isec,isec,isec);
925     if (isec<18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),vecL[isec],dummy);
926     if (isec>=18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),dummy,vecL[isec]);
927     padFitLCE->Add(padFitTmpCE);
928   }
929   //
930   padFitLCE->SetName("CELocal");
931   preprocesor->AddComponent(padFitLCE->Clone());
932   preprocesor->DumpToFile("cetmean.root");   // add it to the file
933   //
934   // write data to array
935   //
936   fFitCELocal  = new TObjArray(6); 
937   fFitCEGlobal = new TObjArray(8); 
938   for (Int_t ipar=0; ipar<8;ipar++){
939     //
940     fFitCEGlobal->AddAt(new TVectorD(36),ipar);
941     TVectorD &fvecG = *((TVectorD*)fFitCEGlobal->At(ipar));
942     for (Int_t isec=0; isec<36;isec++){      
943       if (isec<18)  fvecG[isec]=vecG[0][ipar];
944       if (isec>=18) fvecG[isec]=vecG[1][ipar];    
945     }
946     if (ipar<6){
947       fFitCELocal->AddAt(new TVectorD(36),ipar);
948       TVectorD &fvecL = *((TVectorD*)fFitCELocal->At(ipar));
949       for (Int_t isec=0; isec<36;isec++){      
950         fvecL[isec]=vecL[isec][ipar];
951       }
952     }
953   }
954 }
955
956 void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){
957   //
958   // make a aliases of pad variables
959   //
960   chain->SetAlias("side","(-1+(sector%36<18)*2)");
961   chain->SetAlias("sideA","(sector%36<18)");
962   chain->SetAlias("sideC","(sector%36>=18)");
963   chain->SetAlias("isin","(sector<36)");
964   chain->SetAlias("deltaT","CETmean.fElements-PulserTmean.fElements");
965   chain->SetAlias("timeP","PulserTmean.fElements");
966   chain->SetAlias("Gy","(gy.fElements/500.)");
967   chain->SetAlias("Gx","(gx.fElements/500.)");
968   chain->SetAlias("Lx","(lx.fElements-133)/100.");   // lx in meters
969   chain->SetAlias("Ly","(ly.fElements)/100.");
970   chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)");
971   chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)");
972 }
973
974
975
976
977 void AliTPCkalmanAlign::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
978   //
979   // Update parameters and covariance - with one measurement
980   //
981   const Int_t knMeas=1;
982   Int_t knElem=vecXk.GetNrows();
983  
984   TMatrixD mat1(knElem,knElem);            // update covariance matrix
985   TMatrixD matHk(1,knElem);        // vector to mesurement
986   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
987   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
988   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
989   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
990   TMatrixD covXk2(knElem,knElem);  // helper matrix
991   TMatrixD covXk3(knElem,knElem);  // helper matrix
992   TMatrixD vecZk(1,1);
993   TMatrixD measR(1,1);
994   vecZk(0,0)=delta;
995   measR(0,0)=sigma*sigma;
996   //
997   // reset matHk
998   for (Int_t iel=0;iel<knElem;iel++) 
999     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
1000   //mat1
1001   for (Int_t iel=0;iel<knElem;iel++) {
1002     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
1003     mat1(iel,iel)=1;
1004   }
1005   //
1006   matHk(0, s1)=1;
1007   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
1008   matHkT=matHk.T(); matHk.T();
1009   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
1010   matSk.Invert();
1011   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
1012   vecXk += matKk*vecYk;                    //  updated vector 
1013   covXk2= (mat1-(matKk*matHk));
1014   covXk3 =  covXk2*covXk;          
1015   covXk = covXk3;  
1016   Int_t nrows=covXk3.GetNrows();
1017   
1018   for (Int_t irow=0; irow<nrows; irow++)
1019     for (Int_t icol=0; icol<nrows; icol++){
1020       // rounding problems - make matrix again symteric
1021       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; 
1022     }
1023 }
1024
1025
1026 void   AliTPCkalmanAlign::Update1D(TString &input, TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
1027   //
1028   //  Update Parameters
1029   //  input  - variable name description
1030   //  filter - filter string 
1031   //  param  - parameter vector
1032   //  covar  - covariance
1033   //  mean   - value to set
1034   //  sigma  - sigma of measurement 
1035   TObjArray *array0= input.Tokenize("++");
1036   TObjArray *array1= filter.Tokenize("++");
1037   TMatrixD paramM(param.GetNrows(),1);
1038   for (Int_t i=0; i<array0->GetEntries(); i++){paramM(i,0)=param(i);}
1039
1040   for (Int_t i=0; i<array0->GetEntries(); i++){
1041     Bool_t isOK=kTRUE;
1042     TString str(array0->At(i)->GetName());
1043     for (Int_t j=0; j<array1->GetEntries(); j++){
1044       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1045     }
1046     if (isOK) {
1047       AliTPCkalmanAlign::Update1D(mean, sigma, i, paramM, covar);
1048     }
1049   }
1050   for (Int_t i=0; i<array0->GetEntries(); i++){
1051     param(i)=paramM(i,0);
1052   }
1053 }
1054
1055
1056 TString  AliTPCkalmanAlign::FilterFit(TString &input, TString filter, TVectorD &param, TMatrixD & covar){
1057   //
1058   //
1059   //
1060   TObjArray *array0= input.Tokenize("++");
1061   TObjArray *array1= filter.Tokenize("++");
1062   //TString *presult=new TString("(0");
1063   TString result="(0.0";
1064   for (Int_t i=0; i<array0->GetEntries(); i++){
1065     Bool_t isOK=kTRUE;
1066     TString str(array0->At(i)->GetName());
1067     for (Int_t j=0; j<array1->GetEntries(); j++){
1068       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1069     }
1070     if (isOK) {
1071       result+="+"+str;
1072       result+=Form("*(%f)",param[i+1]);
1073       printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1074     }
1075   }
1076   result+="-0.)";
1077   return result;
1078 }
1079
1080