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