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