]> 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
Adding the class for Cosmic events
[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"
108#include "AliExternalTrackParam.h"
e4042305 109#include "AliTPCTracklet.h"
110#include "TH1D.h"
bb6bc8f6 111#include "TH2F.h"
7eaa723e 112#include "TVectorD.h"
e149f26d 113#include "TTreeStream.h"
7eaa723e 114#include "TFile.h"
6f387311 115#include "TTree.h"
e81dc112 116#include "TF1.h"
8b3c60d8 117#include "TGraphErrors.h"
967eae0d 118#include "AliTPCclusterMI.h"
119#include "AliTPCseed.h"
120#include "AliTracker.h"
121#include "TClonesArray.h"
1d82fc56 122#include "AliExternalComparison.h"
9318a5b4 123
8b3c60d8 124
125#include "TTreeStream.h"
9318a5b4 126#include <iostream>
e4042305 127#include <sstream>
9318a5b4 128using namespace std;
129
6f387311 130AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0;
9318a5b4 131ClassImp(AliTPCcalibAlign)
132
6f387311 133
134
135
136AliTPCcalibAlign* AliTPCcalibAlign::Instance()
137{
138 //
139 // Singleton implementation
140 // Returns an instance of this class, it is created if neccessary
141 //
142 if (fgInstance == 0){
143 fgInstance = new AliTPCcalibAlign();
144 }
145 return fgInstance;
146}
147
148
149
150
9318a5b4 151AliTPCcalibAlign::AliTPCcalibAlign()
bb6bc8f6 152 : AliTPCcalibBase(),
153 fDphiHistArray(72*72),
e4042305 154 fDthetaHistArray(72*72),
155 fDyHistArray(72*72),
156 fDzHistArray(72*72),
bb6bc8f6 157 //
158 fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi
159 fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta
160 fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz
161 fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz
162 fDyZHistArray(72*72), // array of residual histograms y -kYz
163 fDzZHistArray(72*72), // array of residual histograms z -kZz
e4042305 164 fFitterArray12(72*72),
165 fFitterArray9(72*72),
6f387311 166 fFitterArray6(72*72),
167 fMatrixArray12(72*72),
168 fMatrixArray9(72*72),
1d82fc56 169 fMatrixArray6(72*72),
170 fCombinedMatrixArray6(72),
171 fCompTracklet(0), // tracklet comparison
172 fNoField(kFALSE)
9318a5b4 173{
174 //
175 // Constructor
176 //
177 for (Int_t i=0;i<72*72;++i) {
178 fPoints[i]=0;
179 }
180}
181
e149f26d 182AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
183 :AliTPCcalibBase(),
184 fDphiHistArray(72*72),
185 fDthetaHistArray(72*72),
186 fDyHistArray(72*72),
187 fDzHistArray(72*72),
bb6bc8f6 188 fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi
189 fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta
190 fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz
191 fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz
192 fDyZHistArray(72*72), // array of residual histograms y -kYz
6f387311 193 fDzZHistArray(72*72), // array of residual histograms z -kZz //
e149f26d 194 fFitterArray12(72*72),
195 fFitterArray9(72*72),
6f387311 196 fFitterArray6(72*72),
197 fMatrixArray12(72*72),
198 fMatrixArray9(72*72),
1d82fc56 199 fMatrixArray6(72*72),
200 fCombinedMatrixArray6(72),
201 fCompTracklet(0), // tracklet comparison
202 fNoField(kFALSE)
e149f26d 203{
204 //
205 // Constructor
206 //
207 SetName(name);
208 SetTitle(title);
209 for (Int_t i=0;i<72*72;++i) {
210 fPoints[i]=0;
211 }
212}
213
bb6bc8f6 214
215AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
216 :AliTPCcalibBase(align),
217 fDphiHistArray(align.fDphiHistArray),
218 fDthetaHistArray(align.fDthetaHistArray),
219 fDyHistArray(align.fDyHistArray),
220 fDzHistArray(align.fDzHistArray),
221 fDyPhiHistArray(align.fDyPhiHistArray), // array of residual histograms y -kYPhi
222 fDzThetaHistArray(align.fDzThetaHistArray), // array of residual histograms z-z -kZTheta
223 fDphiZHistArray(align.fDphiZHistArray), // array of residual histograms phi -kPhiz
224 fDthetaZHistArray(align.fDthetaZHistArray), // array of residual histograms theta -kThetaz
225 fDyZHistArray(align.fDyZHistArray), // array of residual histograms y -kYz
226 fDzZHistArray(align.fDzZHistArray), // array of residual histograms z -kZz
227 //
228 fFitterArray12(align.fFitterArray12),
229 fFitterArray9(align.fFitterArray9),
6f387311 230 fFitterArray6(align.fFitterArray6),
231
232 fMatrixArray12(align.fMatrixArray12),
233 fMatrixArray9(align.fMatrixArray9),
1d82fc56 234 fMatrixArray6(align.fMatrixArray6),
235 fCombinedMatrixArray6(align.fCombinedMatrixArray6),
236 fCompTracklet(align.fCompTracklet), // tracklet comparison
237 fNoField(align.fNoField)
bb6bc8f6 238{
239 //
240 // copy constructor - copy also the content
241 //
242 TH1 * his = 0;
243 TObjArray * arr0=0;
244 const TObjArray *arr1=0;
245 for (Int_t index =0; index<72*72; index++){
246 for (Int_t iarray=0;iarray<10; iarray++){
247 if (iarray==kY){
248 arr0 = &fDyHistArray;
249 arr1 = &align.fDyHistArray;
250 }
251 if (iarray==kZ){
252 arr0 = &fDzHistArray;
253 arr1 = &align.fDzHistArray;
254 }
255 if (iarray==kPhi){
256 arr0 = &fDphiHistArray;
257 arr1 = &align.fDphiHistArray;
258 }
259 if (iarray==kTheta){
260 arr0 = &fDthetaHistArray;
261 arr1 = &align.fDthetaHistArray;
262 }
263 if (iarray==kYz){
264 arr0 = &fDyZHistArray;
265 arr1 = &align.fDyZHistArray;
266 }
267 if (iarray==kZz){
268 arr0 = &fDzZHistArray;
269 arr1 = &align.fDzZHistArray;
270 }
271 if (iarray==kPhiZ){
272 arr0 = &fDphiZHistArray;
273 arr1 = &align.fDphiZHistArray;
274 }
275 if (iarray==kThetaZ){
276 arr0 = &fDthetaZHistArray;
277 arr1 = &align.fDthetaZHistArray;
278 }
279
280 if (iarray==kYPhi){
281 arr0 = &fDyPhiHistArray;
282 arr1 = &align.fDyPhiHistArray;
283 }
284 if (iarray==kZTheta){
285 arr0 = &fDzThetaHistArray;
286 arr1 = &align.fDzThetaHistArray;
287 }
288
289 if (arr1->At(index)) {
290 his = (TH1*)arr1->At(index)->Clone();
291 his->SetDirectory(0);
292 arr0->AddAt(his,index);
293 }
294 }
295 }
296}
297
298
9318a5b4 299AliTPCcalibAlign::~AliTPCcalibAlign() {
300 //
301 // destructor
302 //
303}
304
e4042305 305void AliTPCcalibAlign::Process(AliTPCseed *seed) {
175d237b 306 //
307 //
308 //
309
e4042305 310 TObjArray tracklets=
311 AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
312 kFALSE,20,2);
175d237b 313 // TObjArray trackletsL=
cbc19295 314// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
315// kFALSE,20,2);
316// TObjArray trackletsQ=
317// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
318// kFALSE,20,2);
319// TObjArray trackletsR=
320// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
321// kFALSE,20,2);
e4042305 322 tracklets.SetOwner();
cbc19295 323 // trackletsL.SetOwner();
324// trackletsQ.SetOwner();
325// trackletsR.SetOwner();
e4042305 326 if (tracklets.GetEntries()==2) {
327 AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
328 AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
329 if (t1->GetSector()>t2->GetSector()) {
330 AliTPCTracklet* tmp=t1;
331 t1=t2;
332 t2=tmp;
333 }
334 AliExternalTrackParam *common1=0,*common2=0;
335 if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
967eae0d 336 ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
e4042305 337 delete common1;
338 delete common2;
339 }
cbc19295 340
e4042305 341}
342
7eaa723e 343void AliTPCcalibAlign::Analyze(){
344 //
345 // Analyze function
346 //
347 EvalFitters();
348}
349
350
351void AliTPCcalibAlign::Terminate(){
352 //
353 // Terminate function
354 // call base terminate + Eval of fitters
355 //
108953e9 356 Info("AliTPCcalibAlign","Terminate");
7eaa723e 357 EvalFitters();
358 AliTPCcalibBase::Terminate();
359}
360
361
362
363
e4042305 364void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
365 const AliExternalTrackParam &tp2,
967eae0d 366 const AliTPCseed * seed,
e4042305 367 Int_t s1,Int_t s2) {
368
972cf6f2 369 //
370 //
371 //
8b3c60d8 372 //
9318a5b4 373 //
374 // Process function to fill fitters
375 //
1d82fc56 376 Double_t t1[10],t2[10];
9318a5b4 377 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
7eaa723e 378 Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 379 x1 =tp1.GetX();
380 y1 =tp1.GetY();
381 z1 =tp1.GetZ();
382 Double_t snp1=tp1.GetSnp();
383 dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
384 Double_t tgl1=tp1.GetTgl();
385 // dz/dx = 1/(cos(theta)*cos(phi))
7eaa723e 386 dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1);
387 x2 =tp2.GetX();
9318a5b4 388 y2 =tp2.GetY();
389 z2 =tp2.GetZ();
390 Double_t snp2=tp2.GetSnp();
391 dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
392 Double_t tgl2=tp2.GetTgl();
7eaa723e 393 dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2);
1d82fc56 394 Int_t accept = AcceptTracklet(tp1,tp2);
7eaa723e 395 //
396 //
397 //
1d82fc56 398 if (fStreamLevel>1 && seed){
7eaa723e 399 TTreeSRedirector *cstream = GetDebugStreamer();
400 if (cstream){
401 static TVectorD vec1(5);
402 static TVectorD vec2(5);
403 vec1.SetElements(t1);
404 vec2.SetElements(t2);
405 AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
406 AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
407 (*cstream)<<"Tracklet"<<
1d82fc56 408 "accept="<<accept<<
108953e9 409 "run="<<fRun<< // run number
410 "event="<<fEvent<< // event number
411 "time="<<fTime<< // time stamp of event
412 "trigger="<<fTrigger<< // trigger
413 "mag="<<fMagF<< // magnetic field
6f387311 414 "isOK="<<accept<< // flag - used for alignment
7eaa723e 415 "tp1.="<<p1<<
416 "tp2.="<<p2<<
417 "v1.="<<&vec1<<
418 "v2.="<<&vec2<<
419 "s1="<<s1<<
420 "s2="<<s2<<
421 "\n";
422 }
423 }
1d82fc56 424 if (accept>0) return;
425 t1[0]-=134.;
426 t2[0]-=134.;
427 t1[5]=0; t2[5]=0;
428 t1[6]=TMath::Sqrt(tp1.GetSigmaY2()+tp2.GetSigmaY2()); t2[6]=t1[6];
429 t1[7]=TMath::Sqrt(tp1.GetSigmaZ2()+tp2.GetSigmaZ2()); t2[7]=t1[7];
430 t1[8]=TMath::Sqrt(tp1.GetSigmaSnp2()+tp2.GetSigmaSnp2()); t2[8]=t1[8];
431 t1[9]=TMath::Sqrt(tp1.GetSigmaTgl2()+tp2.GetSigmaTgl2()); t2[9]=t1[9];
ae0ac7be 432
175d237b 433 if (GetDebugLevel()>50) printf("Process track\n");
6f387311 434 if (GetDebugLevel()>50) printf("Filling track\n");
8b3c60d8 435 //
436 // fill resolution histograms - previous cut included
1d82fc56 437 if (seed) ProcessDiff(tp1,tp2, seed,s1,s2);
8b3c60d8 438 FillHisto(tp1,tp2,s1,s2);
6f387311 439 ProcessAlign(t1,t2,s1,s2);
440}
441
442void AliTPCcalibAlign::ProcessAlign(Double_t * t1,
443 Double_t * t2,
444 Int_t s1,Int_t s2){
445 //
446 // Do intersector alignment
8b3c60d8 447 //
9318a5b4 448 Process12(t1,t2,GetOrMakeFitter12(s1,s2));
7eaa723e 449 Process9(t1,t2,GetOrMakeFitter9(s1,s2));
450 Process6(t1,t2,GetOrMakeFitter6(s1,s2));
e81dc112 451 ++fPoints[GetIndex(s1,s2)];
9318a5b4 452}
453
1d82fc56 454void AliTPCcalibAlign::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){
6f387311 455 //
456 // Process the debug streamer tree
457 // Possible to modify selection criteria
458 // Used with entry list
459 //
460 TTreeSRedirector * cstream = new TTreeSRedirector("aligndump.root");
461
462 AliTPCcalibAlign *align = this;
463 //
464 TVectorD * vec1 = 0;
465 TVectorD * vec2 = 0;
1d82fc56 466 AliExternalTrackParam * tp1 = 0;
467 AliExternalTrackParam * tp2 = 0;
6f387311 468 Int_t s1 = 0;
1d82fc56 469 Int_t s2 = 0;
470 Int_t npoints =0;
6f387311 471 {
472 Int_t entries=chainTracklet->GetEntries();
473 for (Int_t i=0; i< entries; i++){
474 chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
475 chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
476 chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
477 chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
478 chainTracklet->GetBranch("s1")->SetAddress(&s1);
479 chainTracklet->GetBranch("s2")->SetAddress(&s2);
480 chainTracklet->GetEntry(i);
1d82fc56 481 if (!vec1) continue;
482 if (!vec2) continue;
483 if (!tp1) continue;
484 if (!tp2) continue;
485 if (!vec1->GetMatrixArray()) continue;
486 if (!vec2->GetMatrixArray()) continue;
487 // make a local copy
488 AliExternalTrackParam par1(*tp1);
489 AliExternalTrackParam par2(*tp2);
490 TVectorD svec1(*vec1);
491 TVectorD svec2(*vec2);
492 //
6f387311 493 if (s1==s2) continue;
1d82fc56 494 if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i, npoints,s1,s2);
495 AliExternalTrackParam cpar1(par1);
496 AliExternalTrackParam cpar2(par2);
497 Constrain1Pt(cpar1,par2,fNoField);
498 Constrain1Pt(cpar2,par1,fNoField);
499 Bool_t acceptComp = kFALSE;
500 if (comp) acceptComp=comp->AcceptPair(&par1,&par2);
501 if (comp) acceptComp&=comp->AcceptPair(&cpar1,&cpar2);
502 //
503 Int_t reject = align->AcceptTracklet(par1,par2);
504 Int_t rejectC =align->AcceptTracklet(cpar1,cpar2);
505
506 if (1||fStreamLevel>0){
6f387311 507 (*cstream)<<"Tracklet"<<
6f387311 508 "s1="<<s1<<
509 "s2="<<s2<<
1d82fc56 510 "reject="<<reject<<
511 "rejectC="<<rejectC<<
512 "acceptComp="<<acceptComp<<
513 "tp1.="<<&par1<<
514 "tp2.="<<&par2<<
515 "ctp1.="<<&cpar1<<
516 "ctp2.="<<&cpar2<<
517 "v1.="<<&svec1<<
518 "v2.="<<&svec2<<
6f387311 519 "\n";
520 }
1d82fc56 521 //
522 if (fNoField){
523 //
524 //
525 }
526 if (acceptComp) comp->Process(&cpar1,&cpar2);
527 //
528 if (reject>0 || rejectC>0) continue;
529 npoints++;
530 align->ProcessTracklets(cpar1,cpar2,0,s1,s2);
531 align->ProcessTracklets(cpar2,cpar1,0,s2,s1);
6f387311 532 }
533 }
534 delete cstream;
535}
536
537
1d82fc56 538Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
539 const AliExternalTrackParam &p2){
6f387311 540
541 //
542 // Accept pair of tracklets?
543 //
544 /*
545 // resolution cuts
1d82fc56 546 TCut cutS0("sqrt(tp2.fC[0]+tp1.fC[0])<0.2");
547 TCut cutS1("sqrt(tp2.fC[2]+tp1.fC[2])<0.2");
548 TCut cutS2("sqrt(tp2.fC[5]+tp1.fC[5])<0.01");
549 TCut cutS3("sqrt(tp2.fC[9]+tp1.fC[9])<0.01");
550 TCut cutS4("sqrt(tp2.fC[14]+tp1.fC[14])<0.25");
6f387311 551 TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
552 //
553 // parameters matching cuts
554 TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
555 TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
556 TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03");
557 TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03");
1d82fc56 558 TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
559 TCut cutPP4("abs(tp1.fP[4]-tp2.fP[4])/sqrt(tp2.fC[14]+tp1.fC[14])<3");
560 TCut cutP=cutP0+cutP1+cutP2+cutP3+cutP4+cutPP4;
6f387311 561 */
562 //
563 // resolution cuts
1d82fc56 564 Int_t reject=0;
6f387311 565 const Double_t *cp1 = p1.GetCovariance();
566 const Double_t *cp2 = p2.GetCovariance();
1d82fc56 567 if (TMath::Sqrt(cp1[0]+cp2[0])>0.2) reject|=1;;
568 if (TMath::Sqrt(cp1[2]+cp2[2])>0.2) reject|=2;
569 if (TMath::Sqrt(cp1[5]+cp2[5])>0.01) reject|=4;
570 if (TMath::Sqrt(cp1[9]+cp2[9])>0.01) reject|=8;
571 if (TMath::Sqrt(cp1[14]+cp2[14])>0.2) reject|=16;
6f387311 572
573 //parameters difference
574 const Double_t *tp1 = p1.GetParameter();
575 const Double_t *tp2 = p2.GetParameter();
1d82fc56 576 if (TMath::Abs(tp1[0]-tp2[0])>0.6) reject|=32;
577 if (TMath::Abs(tp1[1]-tp2[1])>0.6) reject|=64;
578 if (TMath::Abs(tp1[2]-tp2[2])>0.03) reject|=128;
579 if (TMath::Abs(tp1[3]-tp2[3])>0.03) reject|=526;
580 if (TMath::Abs(tp1[4]-tp2[4])>0.4) reject|=1024;
581 if (TMath::Abs(tp1[4]-tp2[4])/TMath::Sqrt(cp1[14]+cp2[14])>4) reject|=2048;
582
6f387311 583 //
1d82fc56 584 if (TMath::Abs(tp2[1])>235) reject|=2*4096;
585
586 if (fNoField){
587
588 }
589
590 return reject;
6f387311 591}
592
593
594
967eae0d 595void AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
596 const AliExternalTrackParam &t2,
597 const AliTPCseed *seed,
598 Int_t s1,Int_t s2)
599{
600 //
601 // Process local residuals function
602 //
603 TVectorD vecX(160);
604 TVectorD vecY(160);
605 TVectorD vecZ(160);
606 TVectorD vecClY(160);
607 TVectorD vecClZ(160);
608 TClonesArray arrCl("AliTPCclusterMI",160);
609 arrCl.ExpandCreateFast(160);
610 Int_t count1=0, count2=0;
1d82fc56 611
967eae0d 612 for (Int_t i=0;i<160;++i) {
613 AliTPCclusterMI *c=seed->GetClusterPointer(i);
614 vecX[i]=0;
615 vecY[i]=0;
616 vecZ[i]=0;
617 if (!c) continue;
618 AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
619 if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
620 vecClY[i] = c->GetY();
621 vecClZ[i] = c->GetZ();
622 cl=*c;
623 const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
624 if (c->GetDetector()==s1) ++count1;
625 if (c->GetDetector()==s2) ++count2;
626 Double_t gxyz[3],xyz[3];
627 t1.GetXYZ(gxyz);
628 Float_t bz = AliTracker::GetBz(gxyz);
629 par->GetYAt(c->GetX(), bz, xyz[1]);
630 par->GetZAt(c->GetX(), bz, xyz[2]);
631 vecX[i] = c->GetX();
632 vecY[i]= xyz[1];
633 vecZ[i]= xyz[2];
634 }
635 //
636 //
637 if (fStreamLevel>5){
638 //
639 // huge output - cluster residuals to be investigated
640 //
641 TTreeSRedirector *cstream = GetDebugStreamer();
967eae0d 642 AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
643 AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
1c1a1176 644 /*
645
646 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");
647
648 */
649
967eae0d 650 if (cstream){
651 (*cstream)<<"Track"<<
108953e9 652 "run="<<fRun<< // run number
653 "event="<<fEvent<< // event number
654 "time="<<fTime<< // time stamp of event
655 "trigger="<<fTrigger<< // trigger
656 "mag="<<fMagF<< // magnetic field
967eae0d 657 "Cl.="<<&arrCl<<
658 //"tp0.="<<p0<<
659 "tp1.="<<p1<<
660 "tp2.="<<p2<<
661 "vtX.="<<&vecX<<
662 "vtY.="<<&vecY<<
663 "vtZ.="<<&vecZ<<
664 "vcY.="<<&vecClY<<
665 "vcZ.="<<&vecClZ<<
666 "s1="<<s1<<
667 "s2="<<s2<<
668 "c1="<<count1<<
669 "c2="<<count2<<
670 "\n";
671 }
672 }
673}
674
675
676
677
7eaa723e 678void AliTPCcalibAlign::Process12(const Double_t *t1,
679 const Double_t *t2,
9318a5b4 680 TLinearFitter *fitter) {
681 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
682 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
683 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
684 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
685 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
686 //
f8a2dcfb 687 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
688 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
689 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
690
691
692
7eaa723e 693 const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
694 const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 695
1d82fc56 696 //
697 Double_t sy = t1[6];
698 Double_t sz = t1[7];
699 Double_t sdydx = t1[8];
700 Double_t sdzdx = t1[9];
9318a5b4 701
702 Double_t p[12];
703 Double_t value;
704
705 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
706 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
707 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
708 for (Int_t i=0; i<12;i++) p[i]=0.;
709 p[3+0] = x1; // a10
710 p[3+1] = y1; // a11
711 p[3+2] = z1; // a12
712 p[9+1] = 1.; // a13
713 p[0+1] = y1*dydx2; // a01
714 p[0+2] = z1*dydx2; // a02
715 p[9+0] = dydx2; // a03
716 value = y2;
717 fitter->AddPoint(p,value,sy);
718
719 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
720 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
721 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
722 for (Int_t i=0; i<12;i++) p[i]=0.;
723 p[6+0] = x1; // a20
724 p[6+1] = y1; // a21
725 p[6+2] = z1; // a22
726 p[9+2] = 1.; // a23
727 p[0+1] = y1*dzdx2; // a01
728 p[0+2] = z1*dzdx2; // a02
729 p[9+0] = dzdx2; // a03
730 value = z2;
731 fitter->AddPoint(p,value,sz);
732
733 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
734 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
735 for (Int_t i=0; i<12;i++) p[i]=0.;
736 p[3+0] = 1.; // a10
737 p[3+1] = dydx1; // a11
738 p[3+2] = dzdx1; // a12
739 p[0+0] = -dydx2; // a00
740 p[0+1] = -dydx1*dydx2; // a01
741 p[0+2] = -dzdx1*dydx2; // a02
742 value = 0.;
743 fitter->AddPoint(p,value,sdydx);
744
745 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
746 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
747 for (Int_t i=0; i<12;i++) p[i]=0.;
748 p[6+0] = 1; // a20
749 p[6+1] = dydx1; // a21
750 p[6+2] = dzdx1; // a22
751 p[0+0] = -dzdx2; // a00
752 p[0+1] = -dydx1*dzdx2; // a01
753 p[0+2] = -dzdx1*dzdx2; // a02
754 value = 0.;
755 fitter->AddPoint(p,value,sdzdx);
756}
757
758void AliTPCcalibAlign::Process9(Double_t *t1,
759 Double_t *t2,
760 TLinearFitter *fitter) {
761 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
762 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
763 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
764 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
765 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
766 //
f8a2dcfb 767 // a00 a01 a02 a03 1 p[0] p[1] p[6]
768 // a10 a11 a12 a13 ==> p[2] 1 p[3] p[7]
769 // a20 a21 a21 a23 p[4] p[5] 1 p[8]
770
771
9318a5b4 772 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 773 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
1d82fc56 774 //
775 Double_t sy = t1[6];
776 Double_t sz = t1[7];
777 Double_t sdydx = t1[8];
778 Double_t sdzdx = t1[9];
f8a2dcfb 779 //
9318a5b4 780 Double_t p[12];
781 Double_t value;
782
783 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
784 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
f8a2dcfb 785 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
9318a5b4 786 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 787 p[2] += x1; // a10
788 //p[] +=1; // a11
789 p[3] += z1; // a12
790 p[7] += 1; // a13
791 p[0] += y1*dydx2; // a01
792 p[1] += z1*dydx2; // a02
793 p[6] += dydx2; // a03
794 value = y2-y1; //-a11
9318a5b4 795 fitter->AddPoint(p,value,sy);
f8a2dcfb 796 //
9318a5b4 797 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
798 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
f8a2dcfb 799 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
9318a5b4 800 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 801 p[4] += x1; // a20
802 p[5] += y1; // a21
803 //p[] += z1; // a22
804 p[8] += 1.; // a23
805 p[0] += y1*dzdx2; // a01
806 p[1] += z1*dzdx2; // a02
807 p[6] += dzdx2; // a03
808 value = z2-z1; //-a22
9318a5b4 809 fitter->AddPoint(p,value,sz);
810
811 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
812 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
813 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 814 p[2] += 1.; // a10
815 //p[] += dydx1; // a11
816 p[3] += dzdx1; // a12
817 //p[] += -dydx2; // a00
818 p[0] += -dydx1*dydx2; // a01
819 p[1] += -dzdx1*dydx2; // a02
820 value = -dydx1+dydx2; // -a11 + a00
9318a5b4 821 fitter->AddPoint(p,value,sdydx);
822
823 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
824 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
825 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 826 p[4] += 1; // a20
827 p[5] += dydx1; // a21
828 //p[] += dzdx1; // a22
829 //p[] += -dzdx2; // a00
830 p[0] += -dydx1*dzdx2; // a01
831 p[1] += -dzdx1*dzdx2; // a02
832 value = -dzdx1+dzdx2; // -a22 + a00
9318a5b4 833 fitter->AddPoint(p,value,sdzdx);
834}
835
836void AliTPCcalibAlign::Process6(Double_t *t1,
837 Double_t *t2,
838 TLinearFitter *fitter) {
839 // x2 = 1 *x1 +-a01*y1 + 0 +a03
840 // y2 = a01*x1 + 1 *y1 + 0 +a13
841 // z2 = a20*x1 + a21*y1 + 1 *z1 +a23
842 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
843 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
844 //
f8a2dcfb 845 // a00 a01 a02 a03 1 -p[0] 0 p[3]
846 // a10 a11 a12 a13 ==> p[0] 1 0 p[4]
847 // a20 a21 a21 a23 p[1] p[2] 1 p[5]
848
9318a5b4 849 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 850 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 851
1d82fc56 852 //
853 Double_t sy = t1[6];
854 Double_t sz = t1[7];
855 Double_t sdydx = t1[8];
856 Double_t sdzdx = t1[9];
9318a5b4 857
858 Double_t p[12];
859 Double_t value;
f8a2dcfb 860 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
861 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
9318a5b4 862 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
863 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 864 p[0] += x1; // a10
865 //p[] +=1; // a11
866 //p[] += z1; // a12
867 p[4] += 1; // a13
868 p[0] += -y1*dydx2; // a01
869 //p[] += z1*dydx2; // a02
870 p[3] += dydx2; // a03
871 value = y2-y1; //-a11
9318a5b4 872 fitter->AddPoint(p,value,sy);
f8a2dcfb 873 //
874 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
875 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
876 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
9318a5b4 877 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 878 p[1] += x1; // a20
879 p[2] += y1; // a21
880 //p[] += z1; // a22
881 p[5] += 1.; // a23
882 p[0] += -y1*dzdx2; // a01
883 //p[] += z1*dzdx2; // a02
884 p[3] += dzdx2; // a03
885 value = z2-z1; //-a22
9318a5b4 886 fitter->AddPoint(p,value,sz);
887
888 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
889 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
890 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 891 p[0] += 1.; // a10
1d82fc56 892 //p[] += dydx1; // a11
f8a2dcfb 893 //p[] += dzdx1; // a12
894 //p[] += -dydx2; // a00
1d82fc56 895 //p[0] += dydx1*dydx2; // a01 FIXME- 0912 MI
f8a2dcfb 896 //p[] += -dzdx1*dydx2; // a02
897 value = -dydx1+dydx2; // -a11 + a00
9318a5b4 898 fitter->AddPoint(p,value,sdydx);
899
900 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
901 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
902 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 903 p[1] += 1; // a20
1d82fc56 904 // p[2] += dydx1; // a21 FIXME- 0912 MI
f8a2dcfb 905 //p[] += dzdx1; // a22
906 //p[] += -dzdx2; // a00
1d82fc56 907 //p[0] += dydx1*dzdx2; // a01 FIXME- 0912 MI
f8a2dcfb 908 //p[] += -dzdx1*dzdx2; // a02
909 value = -dzdx1+dzdx2; // -a22 + a00
9318a5b4 910 fitter->AddPoint(p,value,sdzdx);
911}
912
7eaa723e 913
914
915
916void AliTPCcalibAlign::EvalFitters() {
917 //
918 // Analyze function
919 //
920 // Perform the fitting using linear fitters
921 //
922 Int_t kMinPoints =50;
9318a5b4 923 TLinearFitter *f;
7eaa723e 924 TFile fff("alignDebug.root","recreate");
9318a5b4 925 for (Int_t s1=0;s1<72;++s1)
7eaa723e 926 for (Int_t s2=0;s2<72;++s2){
e81dc112 927 if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
928 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
0ebabeb6 929 if (f->Eval()!=0) {
9318a5b4 930 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
7eaa723e 931 f->Write(Form("f12_%d_%d",s1,s2));
932 }else{
933 f->Write(Form("f12_%d_%d",s1,s2));
9318a5b4 934 }
935 }
e81dc112 936 if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
937 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
0ebabeb6 938 if (f->Eval()!=0) {
7eaa723e 939 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
940 }else{
941 f->Write(Form("f9_%d_%d",s1,s2));
942 }
943 }
e81dc112 944 if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
945 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
972cf6f2 946 if (f->Eval()!=0) {
7eaa723e 947 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
948 }else{
949 f->Write(Form("f6_%d_%d",s1,s2));
950 }
951 }
952 }
6f387311 953 TMatrixD mat(4,4);
954 for (Int_t s1=0;s1<72;++s1)
955 for (Int_t s2=0;s2<72;++s2){
956 if (GetTransformation12(s1,s2,mat)){
957 fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2));
958 }
959 if (GetTransformation9(s1,s2,mat)){
960 fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2));
961 }
962 if (GetTransformation6(s1,s2,mat)){
963 fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2));
964 }
9318a5b4 965 }
6f387311 966 //this->Write("align");
967
9318a5b4 968}
969
972cf6f2 970TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
971 //
972 // get or make fitter - general linear transformation
973 //
e81dc112 974 static Int_t counter12=0;
975 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 976 TLinearFitter * fitter = GetFitter12(s1,s2);
977 if (fitter) return fitter;
e81dc112 978 // 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]");
979 fitter =new TLinearFitter(&f12,"");
6f387311 980 fitter->StoreData(kFALSE);
e81dc112 981 fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
982 counter12++;
983 if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter12<<endl;
972cf6f2 984 return fitter;
985}
986
987TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
988 //
989 //get or make fitter - general linear transformation - no scaling
990 //
e81dc112 991 static Int_t counter9=0;
992 static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
972cf6f2 993 TLinearFitter * fitter = GetFitter9(s1,s2);
994 if (fitter) return fitter;
e81dc112 995 // fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
996 fitter =new TLinearFitter(&f9,"");
6f387311 997 fitter->StoreData(kFALSE);
972cf6f2 998 fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
e81dc112 999 counter9++;
1000 if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter9<<endl;
972cf6f2 1001 return fitter;
1002}
1003
1004TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
1005 //
1006 // get or make fitter - 6 paramater linear tranformation
1007 // - no scaling
1008 // - rotation x-y
1009 // - tilting x-z, y-z
e81dc112 1010 static Int_t counter6=0;
1011 static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
972cf6f2 1012 TLinearFitter * fitter = GetFitter6(s1,s2);
1013 if (fitter) return fitter;
e81dc112 1014 // fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1015 fitter=new TLinearFitter(&f6,"");
bb6bc8f6 1016 fitter->StoreData(kTRUE);
972cf6f2 1017 fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
e81dc112 1018 counter6++;
1019 if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<" : "<<counter6<<endl;
972cf6f2 1020 return fitter;
1021}
1022
1023
1024
1025
1026
9318a5b4 1027Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
972cf6f2 1028 //
1029 // GetTransformation matrix - 12 paramaters - generael linear transformation
1030 //
9318a5b4 1031 if (!GetFitter12(s1,s2))
1032 return false;
1033 else {
1034 TVectorD p(12);
9318a5b4 1035 GetFitter12(s1,s2)->GetParameters(p);
9318a5b4 1036 a.ResizeTo(4,4);
972cf6f2 1037 a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
1038 a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
1039 a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
1040 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
9318a5b4 1041 return true;
1042 }
1043}
1044
1045Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
972cf6f2 1046 //
1047 // GetTransformation matrix - 9 paramaters - general linear transformation
1048 // No scaling
1049 //
9318a5b4 1050 if (!GetFitter9(s1,s2))
1051 return false;
1052 else {
1053 TVectorD p(9);
1054 GetFitter9(s1,s2)->GetParameters(p);
1055 a.ResizeTo(4,4);
f8a2dcfb 1056 a[0][0]=1; a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
1057 a[1][0]=p[2]; a[1][1]=1; a[1][2]=p[3]; a[1][3]=p[7];
1058 a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1; a[2][3]=p[8];
972cf6f2 1059 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
9318a5b4 1060 return true;
1061 }
1062}
1063
1064Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
972cf6f2 1065 //
1066 // GetTransformation matrix - 6 paramaters
1067 // 3 translation
1068 // 1 rotation -x-y
1069 // 2 tilting x-z y-z
9318a5b4 1070 if (!GetFitter6(s1,s2))
1071 return false;
1072 else {
1073 TVectorD p(6);
9318a5b4 1074 GetFitter6(s1,s2)->GetParameters(p);
9318a5b4 1075 a.ResizeTo(4,4);
f8a2dcfb 1076 a[0][0]=1; a[0][1]=-p[0];a[0][2]=0; a[0][3]=p[3];
1077 a[1][0]=p[0]; a[1][1]=1; a[1][2]=0; a[1][3]=p[4];
1078 a[2][0]=p[1]; a[2][1]=p[2]; a[2][2]=1; a[2][3]=p[5];
972cf6f2 1079 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
9318a5b4 1080 return true;
1081 }
1082}
972cf6f2 1083
1084void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
1085 const AliExternalTrackParam &tp2,
1086 Int_t s1,Int_t s2) {
1087 //
1088 // Fill residual histograms
8b3c60d8 1089 // Innner-Outer
1090 // Left right - x-y
1091 // A-C side
1092 if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0) {
972cf6f2 1093 GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
1094 GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1095 GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
1096 GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
bb6bc8f6 1097 //
1098 GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
1099 GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1100 GetHisto(kYz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetY()-tp2.GetY());
1101 GetHisto(kZz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetZ()-tp2.GetZ());
1102 //
1103 GetHisto(kYPhi,s1,s2,kTRUE)->Fill(tp1.GetSnp(),tp1.GetY()-tp2.GetY());
1104 GetHisto(kZTheta,s1,s2,kTRUE)->Fill(tp1.GetTgl(),tp1.GetZ()-tp2.GetZ());
1105
1106
972cf6f2 1107 }
1108}
1109
1110
1111
1112TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
1113{
1114 //
1115 // return specified residual histogram - it is only QA
1116 // if force specified the histogram and given histogram is not existing
1117 // new histogram is created
1118 //
1119 if (GetIndex(s1,s2)>=72*72) return 0;
1120 TObjArray *histoArray=0;
1121 switch (type) {
1122 case kY:
1123 histoArray = &fDyHistArray; break;
1124 case kZ:
1125 histoArray = &fDzHistArray; break;
1126 case kPhi:
1127 histoArray = &fDphiHistArray; break;
1128 case kTheta:
1129 histoArray = &fDthetaHistArray; break;
bb6bc8f6 1130 case kYPhi:
1131 histoArray = &fDyPhiHistArray; break;
1132 case kZTheta:
1133 histoArray = &fDzThetaHistArray; break;
1134 case kYz:
1135 histoArray = &fDyZHistArray; break;
1136 case kZz:
1137 histoArray = &fDzZHistArray; break;
1138 case kPhiZ:
1139 histoArray = &fDphiZHistArray; break;
1140 case kThetaZ:
1141 histoArray = &fDthetaZHistArray; break;
972cf6f2 1142 }
1143 TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
1144 if (histo) return histo;
1145 if (force==kFALSE) return 0;
1146 //
1147 stringstream name;
1148 stringstream title;
1149 switch (type) {
1150 case kY:
1151 name<<"hist_y_"<<s1<<"_"<<s2;
1152 title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1153 histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1154 break;
1155 case kZ:
1156 name<<"hist_z_"<<s1<<"_"<<s2;
1157 title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1158 histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1159 break;
1160 case kPhi:
1161 name<<"hist_phi_"<<s1<<"_"<<s2;
1162 title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1163 histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1164 break;
1165 case kTheta:
1166 name<<"hist_theta_"<<s1<<"_"<<s2;
1167 title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1168 histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1169 break;
bb6bc8f6 1170 //
1171 //
1172 case kYPhi:
1173 name<<"hist_yphi_"<<s1<<"_"<<s2;
1174 title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1175 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,128,-0.3,0.3); // +/- 3 mm
1176 break;
1177 case kZTheta:
1178 name<<"hist_ztheta_"<<s1<<"_"<<s2;
1179 title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1180 histo = new TH2F(name.str().c_str(),title.str().c_str(),128,20,-1,1,-0.3,0.3); // +/- 3 mm
1181 break;
1182 //
1183 //
1184 //
1185 case kYz:
1186 name<<"hist_yz_"<<s1<<"_"<<s2;
1187 title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1188 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1189 break;
1190 case kZz:
1191 name<<"hist_zz_"<<s1<<"_"<<s2;
1192 title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1193 histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1194 break;
1195 case kPhiZ:
1196 name<<"hist_phiz_"<<s1<<"_"<<s2;
1197 title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1198 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1199 break;
1200 case kThetaZ:
1201 name<<"hist_thetaz_"<<s1<<"_"<<s2;
1202 title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1203 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1204 break;
1205
1206
972cf6f2 1207 }
1208 histo->SetDirectory(0);
1209 histoArray->AddAt(histo,GetIndex(s1,s2));
1210 return histo;
1211}
8b3c60d8 1212
1213TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
1214 Int_t i0, Int_t i1, FitType type)
1215{
1216 //
1217 //
1218 //
1219 TMatrixD mat;
6a77c9f1 1220 //TObjArray *fitArray=0;
8b3c60d8 1221 Double_t xsec[1000];
1222 Double_t ysec[1000];
1223 Int_t npoints=0;
1224 for (Int_t isec = sec0; isec<=sec1; isec++){
1225 Int_t isec2 = (isec+dsec)%72;
1226 switch (type) {
1227 case k6:
1228 GetTransformation6(isec,isec2,mat);break;
1229 case k9:
1230 GetTransformation9(isec,isec2,mat);break;
1231 case k12:
1232 GetTransformation12(isec,isec2,mat);break;
1233 }
1234 xsec[npoints]=isec;
1235 ysec[npoints]=mat(i0,i1);
1236 ++npoints;
1237 }
1238 TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1239 Char_t name[1000];
1240 sprintf(name,"Mat[%d,%d] Type=%d",i0,i1,type);
1241 gr->SetName(name);
1242 return gr;
1243}
1244
1245void AliTPCcalibAlign::MakeTree(const char *fname){
1246 //
1247 // make tree with alignment cosntant -
1248 // For QA visualization
1249 //
ae0ac7be 1250 /*
1251 TFile f("CalibObjects.root");
1252 TObjArray *array = (TObjArray*)f.Get("TPCCalib");
1253 AliTPCcalibAlign *alignTPC = (AliTPCcalibAlign *)array->At(0);
1254 alignTPC->MakeTree("alignTree.root");
1255 TFile falign("alignTree.root");
1256 Align->Draw("dy")
1257 */
8b3c60d8 1258 const Int_t kMinPoints=50;
1259 TTreeSRedirector cstream(fname);
1260 for (Int_t s1=0;s1<72;++s1)
1261 for (Int_t s2=0;s2<72;++s2){
8b3c60d8 1262 TMatrixD m6;
1d82fc56 1263 TMatrixD m6FX;
8b3c60d8 1264 TMatrixD m9;
1265 TMatrixD m12;
8b3c60d8 1266 Double_t dy=0, dz=0, dphi=0,dtheta=0;
1267 Double_t sy=0, sz=0, sphi=0,stheta=0;
1268 Double_t ny=0, nz=0, nphi=0,ntheta=0;
6f387311 1269 Double_t chi2v12=0, chi2v9=0, chi2v6=0;
1270 Int_t npoints=0;
1d82fc56 1271 TLinearFitter * fitter = 0;
6f387311 1272 if (fPoints[GetIndex(s1,s2)]>kMinPoints){
1273 //
1274 //
1275 //
6f387311 1276 fitter = GetFitter12(s1,s2);
1277 npoints = fitter->GetNpoints();
1278 chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1279 //
1280 fitter = GetFitter9(s1,s2);
1281 npoints = fitter->GetNpoints();
1282 chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1283 //
1284 fitter = GetFitter6(s1,s2);
1285 npoints = fitter->GetNpoints();
1286 chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1d82fc56 1287
1288 //
6f387311 1289 GetTransformation6(s1,s2,m6);
1290 GetTransformation9(s1,s2,m9);
1291 GetTransformation12(s1,s2,m12);
1d82fc56 1292 //
1293 fitter = GetFitter6(s1,s2);
1294 fitter->FixParameter(3,0);
1295 fitter->Eval();
1296 GetTransformation6(s1,s2,m6FX);
1297 //
6f387311 1298 TH1 * his=0;
1299 his = GetHisto(kY,s1,s2);
1300 if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1301 his = GetHisto(kZ,s1,s2);
1302 if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1303 his = GetHisto(kPhi,s1,s2);
1304 if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
1305 his = GetHisto(kTheta,s1,s2);
1306 if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
1307 //
1d82fc56 1308
6f387311 1309 }
1310
1311 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1312 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
1313 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
1314 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1315 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1316 //
1317 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
1318 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
1319 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
1320
1321 //
1322 //
1323 // dy:-(134*m6.fElements[4]+m6.fElements[7])
1324 //
1325 // dphi:-(m6.fElements[4])
1326 //
1327 // dz:134*m6.fElements[8]+m6.fElements[11]
1328 //
1329 // dtheta:m6.fElements[8]
8b3c60d8 1330 //
1331 cstream<<"Align"<<
1332 "s1="<<s1<< // reference sector
1333 "s2="<<s2<< // sector to align
1d82fc56 1334 "m6FX.="<<&m6FX<< // tranformation matrix
8b3c60d8 1335 "m6.="<<&m6<< // tranformation matrix
1336 "m9.="<<&m9<< //
1337 "m12.="<<&m12<<
6f387311 1338 "chi2v12="<<chi2v12<<
1339 "chi2v9="<<chi2v9<<
1340 "chi2v6="<<chi2v6<<
967eae0d 1341 // histograms mean RMS and entries
1342 "dy="<<dy<<
8b3c60d8 1343 "sy="<<sy<<
1344 "ny="<<ny<<
1345 "dz="<<dz<<
1346 "sz="<<sz<<
1347 "nz="<<nz<<
1348 "dphi="<<dphi<<
1349 "sphi="<<sphi<<
1350 "nphi="<<nphi<<
1351 "dtheta="<<dtheta<<
1352 "stheta="<<stheta<<
1353 "ntheta="<<ntheta<<
1354 "\n";
1355 }
1356
1357}
ae0ac7be 1358
1359
1360//_____________________________________________________________________
1361Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
1362 //
1363 // merge function
1364 //
1365 if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
1366 if (!list)
1367 return 0;
1368 if (list->IsEmpty())
1369 return 1;
1370
1371 TIterator* iter = list->MakeIterator();
1372 TObject* obj = 0;
1373 iter->Reset();
1374 Int_t count=0;
6f387311 1375 //
1376 TString str1(GetName());
ae0ac7be 1377 while((obj = iter->Next()) != 0)
1378 {
1379 AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
1380 if (entry == 0) continue;
6f387311 1381 if (str1.CompareTo(entry->GetName())!=0) continue;
ae0ac7be 1382 Add(entry);
1383 count++;
1384 }
1385 return count;
1386}
1387
1388
1389void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
1390 //
bb6bc8f6 1391 // Add entry - used for merging of compoents
ae0ac7be 1392 //
ae0ac7be 1393 for (Int_t i=0; i<72;i++){
1394 for (Int_t j=0; j<72;j++){
6f387311 1395 if (align->fPoints[GetIndex(i,j)]<10) continue;
ae0ac7be 1396 fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
ae0ac7be 1397 //
ae0ac7be 1398 //
ae0ac7be 1399 //
bb6bc8f6 1400 for (Int_t itype=0; itype<10; itype++){
1401 TH1 * his0=0, *his1=0;
1402 his0 = GetHisto((HistoType)itype,i,j);
1403 his1 = align->GetHisto((HistoType)itype,i,j);
1404 if (his1){
1405 if (his0) his0->Add(his1);
1406 else {
1407 his0 = GetHisto(kY,i,j,kTRUE);
1408 his0->Add(his1);
1409 }
1410 }
6f387311 1411 }
ae0ac7be 1412 }
1413 }
1414 TLinearFitter *f0=0;
1415 TLinearFitter *f1=0;
1416 for (Int_t i=0; i<72;i++){
6f387311 1417 for (Int_t j=0; j<72;j++){
1418 if (align->fPoints[GetIndex(i,j)]<20) continue;
1419 //
ae0ac7be 1420 //
1421 // fitter12
1422 f0 = GetFitter12(i,j);
bb6bc8f6 1423 f1 = align->GetFitter12(i,j);
6f387311 1424 if (f1 &&f1->Eval()){
1425 if (f0&&f0->GetNpoints()>10) f0->Add(f1);
ae0ac7be 1426 else {
6f387311 1427 fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
ae0ac7be 1428 }
1429 }
1430 //
1431 // fitter9
1432 f0 = GetFitter9(i,j);
bb6bc8f6 1433 f1 = align->GetFitter9(i,j);
6f387311 1434 if (f1&&f1->Eval()){
1435 if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1436 else {
1437 fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
ae0ac7be 1438 }
1439 }
1440 f0 = GetFitter6(i,j);
bb6bc8f6 1441 f1 = align->GetFitter6(i,j);
6f387311 1442 if (f1 &&f1->Eval()){
1443 if (f0&&f0->GetNpoints()>10) f0->Add(f1);
ae0ac7be 1444 else {
6f387311 1445 fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
ae0ac7be 1446 }
1447 }
1448 }
1449 }
1450}
108953e9 1451
6f387311 1452Double_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){
1453 //
1454 // GetTransformed value
1455 //
1456 //
1457 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1458 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
1459 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
1460 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1461 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1462
1463
1464 const TMatrixD * mat = GetTransformation(s1,s2,type);
1465 if (!mat) {
1466 if (value==0) return x1;
1467 if (value==1) return y1;
1468 if (value==2) return z1;
1469 if (value==3) return dydx1;
1470 if (value==4) return dzdx1;
1471 //
1472 if (value==5) return dydx1;
1473 if (value==6) return dzdx1;
1474 return 0;
1475 }
1476 Double_t valT=0;
108953e9 1477
6f387311 1478 if (value==0){
1479 valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
1480 }
1481
1482 if (value==1){
1483 valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
1484 }
1485 if (value==2){
1486 valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
1487 }
1488 if (value==3){
1489 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1490 valT = (*mat)(1,0) +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
1491 valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1492 }
1493
1494 if (value==4){
1495 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1496 valT = (*mat)(2,0) +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
1497 valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1498 }
1499 //
1500 if (value==5){
1501 // onlys shift in angle
1502 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1503 valT = (*mat)(1,0) +(*mat)(1,1)*dydx1;
1504 }
1505
1506 if (value==6){
1507 // only shift in angle
1508 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1509 valT = (*mat)(2,0) +(*mat)(2,1)*dydx1;
1510 }
1511 //
1512 return valT;
1513}
108953e9 1514
1515
1d82fc56 1516void AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){
1517 //
1518 // Update track parameters t1
1519 //
1520 TMatrixD vecXk(5,1); // X vector
1521 TMatrixD covXk(5,5); // X covariance
1522 TMatrixD matHk(1,5); // vector to mesurement
1523 TMatrixD measR(1,1); // measurement error
1524 //TMatrixD matQk(5,5); // prediction noise vector
1525 TMatrixD vecZk(1,1); // measurement
1526 //
1527 TMatrixD vecYk(1,1); // Innovation or measurement residual
1528 TMatrixD matHkT(5,1);
1529 TMatrixD matSk(1,1); // Innovation (or residual) covariance
1530 TMatrixD matKk(5,1); // Optimal Kalman gain
1531 TMatrixD mat1(5,5); // update covariance matrix
1532 TMatrixD covXk2(5,5); //
1533 TMatrixD covOut(5,5);
1534 //
1535 Double_t *param1=(Double_t*) track1.GetParameter();
1536 Double_t *covar1=(Double_t*) track1.GetCovariance();
1537
1538 //
1539 // copy data to the matrix
1540 for (Int_t ipar=0; ipar<5; ipar++){
1541 vecXk(ipar,0) = param1[ipar];
1542 for (Int_t jpar=0; jpar<5; jpar++){
1543 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
1544 }
1545 }
1546 //
1547 //
1548 //
1549 vecZk(0,0) = track2.GetParameter()[4]; // 1/pt measurement from track 2
1550 measR(0,0) = track2.GetCovariance()[14]; // 1/pt measurement error
1551 if (noField) {
1552 measR(0,0)*=0.000000001;
1553 vecZk(0,0)=0.;
1554 }
1555 //
1556 matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0;
1557 matHk(0,3)= 0; matHk(0,4)= 1; // vector to measurement
1558 //
1559 //
1560 //
1561 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
1562 matHkT=matHk.T(); matHk.T();
1563 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
1564 matSk.Invert();
1565 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
1566 vecXk += matKk*vecYk; // updated vector
1567 mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
1568 covXk2 = (mat1-(matKk*matHk));
1569 covOut = covXk2*covXk;
1570 //
1571 //
1572 //
1573 // copy from matrix to parameters
1574 if (0) {
1575 covOut.Print();
1576 vecXk.Print();
1577 covXk.Print();
1578 track1.Print();
1579 track2.Print();
1580 }
1581
1582 for (Int_t ipar=0; ipar<5; ipar++){
1583 param1[ipar]= vecXk(ipar,0) ;
1584 for (Int_t jpar=0; jpar<5; jpar++){
1585 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
1586 }
1587 }
1588
1589}
1590
1591void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){
1592 //
1593 // Global Align -combine the partial alignment of pair of sectors
1594 // minPoints - minimal number of points - don't use sector alignment wit smaller number
1595 // sysError - error added to the alignemnt error
1596 //
1597 AliTPCcalibAlign * align = this;
1598 TMatrixD * arrayAlign[72];
1599 TMatrixD * arrayAlignDiff[72];
1600 //
1601 for (Int_t i=0;i<72; i++) {
1602 TMatrixD * mat = new TMatrixD(4,4);
1603 mat->UnitMatrix();
1604 arrayAlign[i]=mat;
1605 arrayAlignDiff[i]=(TMatrixD*)(mat->Clone());
1606 }
1607
1608 TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root");
1609 for (Int_t iter=0; iter<niter;iter++){
1610 printf("Iter=\t%d\n",iter);
1611 for (Int_t is0=0;is0<72; is0++) {
1612 //
1613 //TMatrixD *mati0 = arrayAlign[is0];
1614 TMatrixD matDiff(4,4);
1615 Double_t sumw=0;
1616 for (Int_t is1=0;is1<72; is1++) {
1617 Bool_t invers=kFALSE;
1618 Int_t npoints=0;
1619 TMatrixD covar;
1620 TVectorD errors;
1621 const TMatrixD *mat = align->GetTransformation(is0,is1,0);
1622 if (mat){
1623 npoints = align->GetFitter6(is0,is1)->GetNpoints();
1624 if (npoints>minPoints){
1625 align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar);
1626 align->GetFitter6(is0,is1)->GetErrors(errors);
1627 }
1628 }
1629 else{
1630 invers=kTRUE;
1631 mat = align->GetTransformation(is1,is0,0);
1632 if (mat) {
1633 npoints = align->GetFitter6(is1,is0)->GetNpoints();
1634 if (npoints>minPoints){
1635 align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar);
1636 align->GetFitter6(is1,is0)->GetErrors(errors);
1637 }
1638 }
1639 }
1640 if (!mat) continue;
1641 if (npoints<minPoints) continue;
1642 //
1643 Double_t weight=1;
1644 if (is1/36>is0/36) weight*=2./3.; //IROC-OROC
1645 if (is1/36<is0/36) weight*=1./3.; //OROC-IROC
1646 if (is1/36==is0/36) weight*=1/3.; //OROC-OROC
1647 if (is1%36!=is0%36) weight*=1/2.; //Not up-down
1648 weight/=(errors[4]*errors[4]+sysError*sysError);
1649 //
1650 //
1651 TMatrixD matT = *mat;
1652 if (invers) matT.Invert();
1653 TMatrixD diffMat= (*(arrayAlign[is1]))*matT;
1654 diffMat-=(*arrayAlign[is0]);
1655 matDiff+=weight*diffMat;
1656 sumw+=weight;
1657
1658 (*cstream)<<"LAlign"<<
1659 "iter="<<iter<<
1660 "s0="<<is0<<
1661 "s1="<<is1<<
1662 "npoints="<<npoints<<
1663 "m60.="<<arrayAlign[is0]<<
1664 "m61.="<<arrayAlign[is1]<<
1665 "m01.="<<&matT<<
1666 "diff.="<<&diffMat<<
1667 "cov.="<<&covar<<
1668 "err.="<<&errors<<
1669 "\n";
1670 }
1671 if (sumw>0){
1672 matDiff*=1/sumw;
1673 matDiff(0,0)=0;
1674 matDiff(1,1)=0;
1675 matDiff(1,1)=0;
1676 matDiff(1,1)=0;
1677 (*arrayAlignDiff[is0]) = matDiff;
1678 }
1679 }
1680 for (Int_t is0=0;is0<72; is0++) {
1681 if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]);
1682 if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]);
1683 //
1684 (*cstream)<<"GAlign"<<
1685 "iter="<<iter<<
1686 "s0="<<is0<<
1687 "m6.="<<arrayAlign[is0]<<
1688 "\n";
1689 }
1690 }
1691 delete cstream;
1692 for (Int_t isec=0;isec<72;isec++){
1693 fCombinedMatrixArray6.AddAt(arrayAlign[isec],isec);
1694 delete arrayAlignDiff[isec];
1695 }
1696}
1697
108953e9 1698
1699/*
1700
1701
1702gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1703gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1704AliXRDPROOFtoolkit tool;
1705TChain * chainTr = tool.MakeChain("align.txt","Track",0,10200);
1706chainTr->Lookup();
6f387311 1707TChain * chainTracklet = tool.MakeChain("align.txt","Tracklet",0,20);
1708chainTracklet->Lookup();
1709
1710
1711TCut cutS0("sqrt(tp2.fC[0]+tp2.fC[0])<0.6");
1712TCut cutS1("sqrt(tp2.fC[2]+tp2.fC[2])<0.6");
1713TCut cutS2("sqrt(tp2.fC[5]+tp2.fC[5])<0.04");
1714TCut cutS3("sqrt(tp2.fC[9]+tp2.fC[9])<0.02");
1715TCut cutS4("sqrt(tp2.fC[14]+tp2.fC[14])<0.5");
1716// resolution cuts
1717TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
1718
1719TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1720TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1721TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.02");
1722TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.01");
1723TCut cutE("abs(tp2.fP[1])<235");
1724TCut cutP=cutP0+cutP1+cutP2+cutP3+cutE;
1725
108953e9 1726
bb6bc8f6 1727//
1728//
6f387311 1729TCut cutA = cutP+cutS;
1730chainTr->Draw(">>listEL",cutA,"entryList");
bb6bc8f6 1731TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1732chainTr->SetEntryList(elist);
1733
1734
1735
1736
108953e9 1737
1738
1739TCut cutS("s1%36==s2%36");
1740
1741TCut cutN("c1>32&&c2>60");
1742TCut cutC0("sqrt(tp2.fC[0])<1");
1743
108953e9 1744TCut cutX("abs(tp2.fX-133.6)<2");
1745
1746TCut cutA = cutP+cutN;
1747
1748
1749TCut cutY("abs(vcY.fElements-vtY.fElements)<0.3&&vcY.fElements!=0")
1750TCut cutZ("abs(vcZ.fElements-vtZ.fElements)<0.3&&vcZ.fElements!=0")
1751
6f387311 1752
1753
1754
1755
1756
1757TCut cutA = cutP+cutS;
1758chainTracklet->Draw(">>listEL",cutA,"entryList");
1759TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1760chainTracklet->SetEntryList(elist);
1761
1762//
1763TVectorD * vec1 = 0;
1764TVectorD * vec2 = 0;
1765AliExternalTrackParam * tp1 = 0;
1766AliExternalTrackParam * tp2 = 0;
1767
1768Int_t s1 = 0;
1769Int_t s2 = 0;
1770chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1771chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1772chainTracklet->GetBranch("s1")->SetAddress(&s1);
1773chainTracklet->GetBranch("s2")->SetAddress(&s2);
1774
1775
1776AliTPCcalibAlign align;
1777{
1778for (Int_t i=0; i< elist->GetN(); i++){
1779//for (Int_t i=0; i<100000; i++){
1780chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
1781chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
1782chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1783chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1784chainTracklet->GetBranch("s1")->SetAddress(&s1);
1785chainTracklet->GetBranch("s2")->SetAddress(&s2);
1786
1787chainTracklet->GetEntry(i);
1788if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i,tentry, s1,s2);
1789//vec1.Print();
1790TLinearFitter * fitter = align.GetOrMakeFitter6(s1,s2);
1791if (fitter) align.Process6(vec1->GetMatrixArray(),vec2->GetMatrixArray(),fitter);
1792}
1793}
1794
108953e9 1795*/