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