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