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