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