]>
Commit | Line | Data |
---|---|---|
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 | |
79 | AliTPCkalmanAlign::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 | ||
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 | |
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 | ||
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 ¶m, 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; | |
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 | ||
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 | } | |
125d3a38 | 249 | } |
250 | ||
251 | ||
252 | ||
253 | ||
254 | void 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.="<<¶mCE[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.="<<¶mCE[is0%36]<< // default CE parameters |
82628455 | 413 | "pceOut0.="<<¶mCE[is0%36+36]<< |
414 | "pceIn1.="<<¶mCE[is1%36]<< | |
415 | "pceOut1.="<<¶mCE[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 | ||
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(); | |
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 | ||
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); | |
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 | ||
755 | void 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 | ||
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 | // | |
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)); | |
922 | // | |
7d85e147 | 923 | 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 | 924 | if (isec<18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),vecL[isec],dummy); |
925 | if (isec>=18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),dummy,vecL[isec]); | |
926 | padFitLCE->Add(padFitTmpCE); | |
927 | } | |
928 | // | |
929 | padFitLCE->SetName("CELocal"); | |
930 | preprocesor->AddComponent(padFitLCE->Clone()); | |
931 | preprocesor->DumpToFile("cetmean.root"); // add it to the file | |
932 | // | |
933 | // write data to array | |
934 | // | |
4486a91f | 935 | fFitCELocal = new TObjArray(6); |
0b736a46 | 936 | fFitCEGlobal = new TObjArray(8); |
937 | for (Int_t ipar=0; ipar<8;ipar++){ | |
4486a91f | 938 | // |
0b736a46 | 939 | fFitCEGlobal->AddAt(new TVectorD(36),ipar); |
4486a91f | 940 | TVectorD &fvecG = *((TVectorD*)fFitCEGlobal->At(ipar)); |
4486a91f | 941 | for (Int_t isec=0; isec<36;isec++){ |
0b736a46 | 942 | if (isec<18) fvecG[isec]=vecG[0][ipar]; |
943 | if (isec>=18) fvecG[isec]=vecG[1][ipar]; | |
944 | } | |
945 | if (ipar<6){ | |
946 | fFitCELocal->AddAt(new TVectorD(36),ipar); | |
947 | TVectorD &fvecL = *((TVectorD*)fFitCELocal->At(ipar)); | |
948 | for (Int_t isec=0; isec<36;isec++){ | |
949 | fvecL[isec]=vecL[isec][ipar]; | |
4486a91f | 950 | } |
951 | } | |
952 | } | |
4486a91f | 953 | } |
125d3a38 | 954 | |
4486a91f | 955 | void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){ |
956 | // | |
957 | // make a aliases of pad variables | |
958 | // | |
959 | chain->SetAlias("side","(-1+(sector%36<18)*2)"); | |
960 | chain->SetAlias("sideA","(sector%36<18)"); | |
961 | chain->SetAlias("sideC","(sector%36>=18)"); | |
962 | chain->SetAlias("isin","(sector<36)"); | |
963 | chain->SetAlias("deltaT","CETmean.fElements-PulserTmean.fElements"); | |
964 | chain->SetAlias("timeP","PulserTmean.fElements"); | |
965 | chain->SetAlias("Gy","(gy.fElements/500.)"); | |
966 | chain->SetAlias("Gx","(gx.fElements/500.)"); | |
967 | chain->SetAlias("Lx","(lx.fElements-133)/100."); // lx in meters | |
968 | chain->SetAlias("Ly","(ly.fElements)/100."); | |
969 | chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)"); | |
970 | chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)"); | |
125d3a38 | 971 | } |