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