Fixes for building of DA (Anshul)
[u/mrichter/AliRoot.git] / TPC / AliTPCkalmanAlign.cxx
CommitLineData
125d3a38 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 //
4486a91f 34 Int_t run=117220;
125d3a38 35 .x ConfigCalibTrain.C(run)
4486a91f 36 calibDB = AliTPCcalibDB::Instance()
125d3a38 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"
4486a91f 52#include "TCut.h"
53#include "TChain.h"
54#include "AliXRDPROOFtoolkit.h"
55#include "TLegend.h"
56#include "TCanvas.h"
125d3a38 57
82628455 58#include "AliTPCcalibDB.h"
59#include "AliTPCCalROC.h"
125d3a38 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"
82628455 67#include "AliTracker.h"
125d3a38 68#include "TFile.h"
69#include "TLinearFitter.h"
70#include "AliTPCcalibAlign.h"
71#include "TH1.h"
72#include "AliTPCCalPad.h"
73#include "AliTPCkalmanAlign.h"
4486a91f 74#include "TStatToolkit.h"
75#include "AliTPCPreprocessorOnline.h"
76#include "TPostScript.h"
0b736a46 77#include "AliGRPObject.h"
125d3a38 78
79AliTPCkalmanAlign::AliTPCkalmanAlign():
80 TNamed(),
81 fCalibAlign(0), // kalman alignnmnt
82 fOriginalAlign(0), // original alignment 0 read for the OCDB
4486a91f 83 fNewAlign(0),
84 fPadTime0(0),
85 fFitCEGlobal(0),
86 fFitCELocal(0)
125d3a38 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
97AliTPCkalmanAlign::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
4486a91f 101 fNewAlign(0),
102 fPadTime0(0),
103 fFitCEGlobal(0),
104 fFitCELocal(0)
125d3a38 105{
106 //
107 // Default constructor
108 //
109 for (Int_t i=0; i<4; i++){
110 fDelta1D[i]=0;
111 fCovar1D[i]=0;
112 }
4486a91f 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 }
125d3a38 119}
120
121void 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
142void 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
156void 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;
5647625c 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
204void 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 }
125d3a38 249}
250
251
252
253
254void AliTPCkalmanAlign::MakeGlobalAlign(){
255 //
256 // Combine all pairs of fitters and make global alignemnt
257 //
0b736a46 258
125d3a38 259 AliTPCkalmanAlign &kalmanAlign=*this;
260 TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
0b736a46 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
82628455 268 //
269 // get ce info
270 //
271 AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
272 TVectorD paramCE[72];
273 TMatrixD covarCE[72];
4486a91f 274 Int_t statUpDown=0; // statistic up down
275 Int_t statLeftRight=0; // statistic left-right
82628455 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
125d3a38 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);
82628455 297 TVectorD stat(5);
125d3a38 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();
82628455 307 rms[0]=his->GetRMS();
308 stat[0]=his->GetEntries();
125d3a38 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();
82628455 314 rms[1]=his->GetRMS();
315 stat[1]=his->GetEntries();
125d3a38 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();
82628455 321 rms[2]=his->GetRMS();
322 stat[2]=his->GetEntries();
125d3a38 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();
82628455 328 rms[3]=his->GetRMS();
329 stat[3]=his->GetEntries();
125d3a38 330 kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
4486a91f 331 if (is1==is0+36) statUpDown+=Int_t(stat[0]);
332 if (is1==is0+35) statLeftRight+=Int_t(stat[0]);
125d3a38 333 }
82628455 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();
4486a91f 344 statUpDown/=36;
345 statLeftRight/=36;
82628455 346 MakeNewAlignment(kTRUE);
0b736a46 347 FitCE();
125d3a38 348 for (Int_t is0=0;is0<72;is0++)
349 for (Int_t is1=0;is1<72;is1++){
82628455 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;
125d3a38 355 //
356 //
357 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
82628455 358 if (his){
359 delta[0]=his->GetMean();
360 rms[0]=his->GetRMS();
361 stat[0]=his->GetEntries();
362 }
125d3a38 363 //
364 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
82628455 365 if (his) {
366 delta[1]=his->GetMean();
367 rms[1]=his->GetRMS();
368 stat[1]=his->GetEntries();
369 }
125d3a38 370 //
371 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
82628455 372 if (his){
373 delta[2] = his->GetMean();
374 rms[2]=his->GetRMS();
375 stat[2]=his->GetEntries();
376 }
125d3a38 377 //
378 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
82628455 379 if (his){
380 delta[3] = his->GetMean();
381 rms[3]=his->GetRMS();
382 stat[3]=his->GetEntries();
383 }
0b736a46 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
125d3a38 397 (*pcstream)<<"kalmanAlignDebug"<<
0b736a46 398 "run="<<run<< // runs
399 "time="<<time<< // time
400 "timeE="<<timeE<< // sart of tun time
401 "timeS="<<timeS<< // end od run time
82628455 402 "bz="<<bz<<
125d3a38 403 "is0="<<is0<<
404 "is1="<<is1<<
0b736a46 405 "delta.="<<&delta<< // alignment deltas
406 "rms.="<<&rms<< // rms
82628455 407 "stat.="<<&stat<<
125d3a38 408 "vec0.="<<&vec[0]<<
409 "vec1.="<<&vec[1]<<
410 "vec2.="<<&vec[2]<<
411 "vec3.="<<&vec[3]<<
0b736a46 412 "pceIn0.="<<&paramCE[is0%36]<< // default CE parameters
82628455 413 "pceOut0.="<<&paramCE[is0%36+36]<<
414 "pceIn1.="<<&paramCE[is1%36]<<
415 "pceOut1.="<<&paramCE[is1%36+36]<<
0b736a46 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
125d3a38 432 "\n";
433 }
4486a91f 434
4486a91f 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
125d3a38 442 delete pcstream;
443}
444
445
446
447
448
449
450void 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
470void 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();
4486a91f 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;
125d3a38 658}
659
660
661
662void 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
697void 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);
82628455 719 //AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
125d3a38 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 }
82628455 734 localTransNew[1]=localTransNew[1]-((*fDelta1D[0])(sector,0));
735 localRot[0] =localRot[0]-(*fDelta1D[2])(sector,0);
125d3a38 736 //
82628455 737 if (pcstream) (*pcstream)<<"alignParams"<<
125d3a38 738 //"idLayer="<<idLayer<<
739 "idModule="<<idModule<<
740 "sector="<<sector<<
741 "olT.="<<&localTrans<<
125d3a38 742 "olR.="<<&localRot<<
82628455 743 "ogT.="<<&localTrans<<
125d3a38 744 "ogR.="<<&globalRot<<
82628455 745 "nlT.="<<&localTransNew<<
746 "nlR.="<<&localRotNew<<
747 "ngT.="<<&localTransNew<<
748 "ngR.="<<&globalRotNew<<
125d3a38 749 "\n";
750 }
4486a91f 751}
752
753
754
755void AliTPCkalmanAlign::DrawAlignmentTrends(){
756 //
757 // Draw trends of alingment variables
758 //
759 /*
0b736a46 760 //1. Create a list of align data
761 //
762 //2. Filter list
763 AliXRDPROOFtoolkit::FilterList("align.list","AliTPCkalmanAlign.root kalmanAlignDebug",0);
764
4486a91f 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};
125d3a38 781
4486a91f 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
807void 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 //
0b736a46 865 fstringG+="Gy++"; // 1 - global y
866 fstringG+="Gx++"; // 2 - global x
4486a91f 867 //
0b736a46 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
4486a91f 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
0b736a46 912 fstringL+="Ly++"; // 4. common slope
913 fstringL+="Ly*isin++"; // 5. delta slope
4486a91f 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++){
7d85e147 919 TCut cutSector=Form("(sector%%36)==%d",isec);
4486a91f 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));
472f0066 922 delete strFitLX;
4486a91f 923 //
7d85e147 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);
4486a91f 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 //
4486a91f 936 fFitCELocal = new TObjArray(6);
0b736a46 937 fFitCEGlobal = new TObjArray(8);
938 for (Int_t ipar=0; ipar<8;ipar++){
4486a91f 939 //
0b736a46 940 fFitCEGlobal->AddAt(new TVectorD(36),ipar);
4486a91f 941 TVectorD &fvecG = *((TVectorD*)fFitCEGlobal->At(ipar));
4486a91f 942 for (Int_t isec=0; isec<36;isec++){
0b736a46 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];
4486a91f 951 }
952 }
953 }
4486a91f 954}
125d3a38 955
4486a91f 956void 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)");
125d3a38 972}
46e89793 973
974
975
976
977void 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
1026void 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
1056TString 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