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