]> 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 alignment macro (Marian)
[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");
175d237b 99
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"
e81dc112 115#include "TF1.h"
8b3c60d8 116#include "TGraphErrors.h"
967eae0d 117#include "AliTPCclusterMI.h"
118#include "AliTPCseed.h"
119#include "AliTracker.h"
120#include "TClonesArray.h"
9318a5b4 121
8b3c60d8 122
123#include "TTreeStream.h"
9318a5b4 124#include <iostream>
e4042305 125#include <sstream>
9318a5b4 126using namespace std;
127
128ClassImp(AliTPCcalibAlign)
129
130AliTPCcalibAlign::AliTPCcalibAlign()
bb6bc8f6 131 : AliTPCcalibBase(),
132 fDphiHistArray(72*72),
e4042305 133 fDthetaHistArray(72*72),
134 fDyHistArray(72*72),
135 fDzHistArray(72*72),
bb6bc8f6 136 //
137 fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi
138 fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta
139 fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz
140 fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz
141 fDyZHistArray(72*72), // array of residual histograms y -kYz
142 fDzZHistArray(72*72), // array of residual histograms z -kZz
e4042305 143 fFitterArray12(72*72),
144 fFitterArray9(72*72),
145 fFitterArray6(72*72)
9318a5b4 146{
147 //
148 // Constructor
149 //
150 for (Int_t i=0;i<72*72;++i) {
151 fPoints[i]=0;
152 }
153}
154
e149f26d 155AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
156 :AliTPCcalibBase(),
157 fDphiHistArray(72*72),
158 fDthetaHistArray(72*72),
159 fDyHistArray(72*72),
160 fDzHistArray(72*72),
bb6bc8f6 161 fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi
162 fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta
163 fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz
164 fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz
165 fDyZHistArray(72*72), // array of residual histograms y -kYz
166 fDzZHistArray(72*72), // array of residual histograms z -kZz
167
e149f26d 168 fFitterArray12(72*72),
169 fFitterArray9(72*72),
170 fFitterArray6(72*72)
171{
172 //
173 // Constructor
174 //
175 SetName(name);
176 SetTitle(title);
177 for (Int_t i=0;i<72*72;++i) {
178 fPoints[i]=0;
179 }
180}
181
bb6bc8f6 182
183AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
184 :AliTPCcalibBase(align),
185 fDphiHistArray(align.fDphiHistArray),
186 fDthetaHistArray(align.fDthetaHistArray),
187 fDyHistArray(align.fDyHistArray),
188 fDzHistArray(align.fDzHistArray),
189 fDyPhiHistArray(align.fDyPhiHistArray), // array of residual histograms y -kYPhi
190 fDzThetaHistArray(align.fDzThetaHistArray), // array of residual histograms z-z -kZTheta
191 fDphiZHistArray(align.fDphiZHistArray), // array of residual histograms phi -kPhiz
192 fDthetaZHistArray(align.fDthetaZHistArray), // array of residual histograms theta -kThetaz
193 fDyZHistArray(align.fDyZHistArray), // array of residual histograms y -kYz
194 fDzZHistArray(align.fDzZHistArray), // array of residual histograms z -kZz
195 //
196 fFitterArray12(align.fFitterArray12),
197 fFitterArray9(align.fFitterArray9),
198 fFitterArray6(align.fFitterArray6)
199
200{
201 //
202 // copy constructor - copy also the content
203 //
204 TH1 * his = 0;
205 TObjArray * arr0=0;
206 const TObjArray *arr1=0;
207 for (Int_t index =0; index<72*72; index++){
208 for (Int_t iarray=0;iarray<10; iarray++){
209 if (iarray==kY){
210 arr0 = &fDyHistArray;
211 arr1 = &align.fDyHistArray;
212 }
213 if (iarray==kZ){
214 arr0 = &fDzHistArray;
215 arr1 = &align.fDzHistArray;
216 }
217 if (iarray==kPhi){
218 arr0 = &fDphiHistArray;
219 arr1 = &align.fDphiHistArray;
220 }
221 if (iarray==kTheta){
222 arr0 = &fDthetaHistArray;
223 arr1 = &align.fDthetaHistArray;
224 }
225 if (iarray==kYz){
226 arr0 = &fDyZHistArray;
227 arr1 = &align.fDyZHistArray;
228 }
229 if (iarray==kZz){
230 arr0 = &fDzZHistArray;
231 arr1 = &align.fDzZHistArray;
232 }
233 if (iarray==kPhiZ){
234 arr0 = &fDphiZHistArray;
235 arr1 = &align.fDphiZHistArray;
236 }
237 if (iarray==kThetaZ){
238 arr0 = &fDthetaZHistArray;
239 arr1 = &align.fDthetaZHistArray;
240 }
241
242 if (iarray==kYPhi){
243 arr0 = &fDyPhiHistArray;
244 arr1 = &align.fDyPhiHistArray;
245 }
246 if (iarray==kZTheta){
247 arr0 = &fDzThetaHistArray;
248 arr1 = &align.fDzThetaHistArray;
249 }
250
251 if (arr1->At(index)) {
252 his = (TH1*)arr1->At(index)->Clone();
253 his->SetDirectory(0);
254 arr0->AddAt(his,index);
255 }
256 }
257 }
258}
259
260
9318a5b4 261AliTPCcalibAlign::~AliTPCcalibAlign() {
262 //
263 // destructor
264 //
265}
266
e4042305 267void AliTPCcalibAlign::Process(AliTPCseed *seed) {
175d237b 268 //
269 //
270 //
271
e4042305 272 TObjArray tracklets=
273 AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
274 kFALSE,20,2);
175d237b 275 // TObjArray trackletsL=
cbc19295 276// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
277// kFALSE,20,2);
278// TObjArray trackletsQ=
279// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
280// kFALSE,20,2);
281// TObjArray trackletsR=
282// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
283// kFALSE,20,2);
e4042305 284 tracklets.SetOwner();
cbc19295 285 // trackletsL.SetOwner();
286// trackletsQ.SetOwner();
287// trackletsR.SetOwner();
e4042305 288 if (tracklets.GetEntries()==2) {
289 AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
290 AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
291 if (t1->GetSector()>t2->GetSector()) {
292 AliTPCTracklet* tmp=t1;
293 t1=t2;
294 t2=tmp;
295 }
296 AliExternalTrackParam *common1=0,*common2=0;
297 if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
967eae0d 298 ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
e4042305 299 delete common1;
300 delete common2;
301 }
cbc19295 302
e4042305 303}
304
7eaa723e 305void AliTPCcalibAlign::Analyze(){
306 //
307 // Analyze function
308 //
309 EvalFitters();
310}
311
312
313void AliTPCcalibAlign::Terminate(){
314 //
315 // Terminate function
316 // call base terminate + Eval of fitters
317 //
108953e9 318 Info("AliTPCcalibAlign","Terminate");
7eaa723e 319 EvalFitters();
320 AliTPCcalibBase::Terminate();
321}
322
323
324
325
e4042305 326void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
327 const AliExternalTrackParam &tp2,
967eae0d 328 const AliTPCseed * seed,
e4042305 329 Int_t s1,Int_t s2) {
330
972cf6f2 331 //
332 //
333 //
8b3c60d8 334 //
9318a5b4 335 //
336 // Process function to fill fitters
337 //
338 Double_t t1[5],t2[5];
339 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
7eaa723e 340 Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 341 x1 =tp1.GetX();
342 y1 =tp1.GetY();
343 z1 =tp1.GetZ();
344 Double_t snp1=tp1.GetSnp();
345 dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
346 Double_t tgl1=tp1.GetTgl();
347 // dz/dx = 1/(cos(theta)*cos(phi))
7eaa723e 348 dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1);
349 x2 =tp2.GetX();
9318a5b4 350 y2 =tp2.GetY();
351 z2 =tp2.GetZ();
352 Double_t snp2=tp2.GetSnp();
353 dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
354 Double_t tgl2=tp2.GetTgl();
7eaa723e 355 dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2);
7eaa723e 356 //
357 //
358 //
359 if (fStreamLevel>1){
360 TTreeSRedirector *cstream = GetDebugStreamer();
361 if (cstream){
362 static TVectorD vec1(5);
363 static TVectorD vec2(5);
364 vec1.SetElements(t1);
365 vec2.SetElements(t2);
366 AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
367 AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
368 (*cstream)<<"Tracklet"<<
108953e9 369 "run="<<fRun<< // run number
370 "event="<<fEvent<< // event number
371 "time="<<fTime<< // time stamp of event
372 "trigger="<<fTrigger<< // trigger
373 "mag="<<fMagF<< // magnetic field
7eaa723e 374 "tp1.="<<p1<<
375 "tp2.="<<p2<<
376 "v1.="<<&vec1<<
377 "v2.="<<&vec2<<
378 "s1="<<s1<<
379 "s2="<<s2<<
380 "\n";
381 }
382 }
8b3c60d8 383 //
384 // Aplly cut selection
385 /*
ae0ac7be 386 TChain * chainalign = tool.MakeChain("align.txt","Tracklet",0,1000000)
387 chainalign->Lookup();
388
8b3c60d8 389 // Cuts to be justified with the debug streamer
390 //
7b18d067 391 TCut c1pt("abs((tp1.fP[4]+tp2.fP[4])*0.5)<3"); // pt cut - OK
392 TCut cdy("abs(tp1.fP[0]-tp2.fP[0])<2");
393 TCut cdz("abs(tp1.fP[1]-tp2.fP[1])<2");
394 TCut cdphi("abs(tp1.fP[2]-tp2.fP[2])<0.02");
395 TCut cdt("abs(tp1.fP[3]-tp2.fP[3])<0.02");
396 TCut cd1pt("abs(tp1.fP[4]-tp2.fP[4])<0.3"); // delta 1/pt cut -OK
8b3c60d8 397 //
398 //
7b18d067 399 TCut acut = c1pt+cdy+cdz+cdphi+cdt+cd1pt;
bb6bc8f6 400 chainalign->Draw(">>listEL",acut,"entryList");
401 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
402 chainalign->SetEntryList(elist);
ae0ac7be 403
8b3c60d8 404 */
405 // 1. pt cut
7b18d067 406 // 2. dy
407 // 3. dz
408 // 4. dphi
409 // 5. dtheta
410 // 6. d1pt
175d237b 411 if (GetDebugLevel()>50) printf("Process track\n");
7b18d067 412 if (TMath::Abs(tp1.GetParameter()[0]-tp2.GetParameter()[0])>2) return;
413 if (TMath::Abs(tp1.GetParameter()[1]-tp2.GetParameter()[1])>2) return;
414 if (TMath::Abs(tp1.GetParameter()[2]-tp2.GetParameter()[2])>0.02) return;
415 if (TMath::Abs(tp1.GetParameter()[3]-tp2.GetParameter()[3])>0.02) return;
416 if (TMath::Abs(tp1.GetParameter()[4]-tp2.GetParameter()[4])>0.3) return;
ae0ac7be 417 if (TMath::Abs((tp1.GetParameter()[4]+tp2.GetParameter()[4])*0.5)>3) return;
418 if (TMath::Abs((tp1.GetParameter()[0]-tp2.GetParameter()[0]))<0.000000001) return;
175d237b 419 if (GetDebugLevel()>50) printf("Filling track\n");
8b3c60d8 420 //
421 // fill resolution histograms - previous cut included
422 FillHisto(tp1,tp2,s1,s2);
423 //
9318a5b4 424 Process12(t1,t2,GetOrMakeFitter12(s1,s2));
7eaa723e 425 Process9(t1,t2,GetOrMakeFitter9(s1,s2));
426 Process6(t1,t2,GetOrMakeFitter6(s1,s2));
967eae0d 427 ProcessDiff(tp1,tp2, seed,s1,s2);
e81dc112 428 ++fPoints[GetIndex(s1,s2)];
9318a5b4 429}
430
967eae0d 431void AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
432 const AliExternalTrackParam &t2,
433 const AliTPCseed *seed,
434 Int_t s1,Int_t s2)
435{
436 //
437 // Process local residuals function
438 //
439 TVectorD vecX(160);
440 TVectorD vecY(160);
441 TVectorD vecZ(160);
442 TVectorD vecClY(160);
443 TVectorD vecClZ(160);
444 TClonesArray arrCl("AliTPCclusterMI",160);
445 arrCl.ExpandCreateFast(160);
446 Int_t count1=0, count2=0;
447 for (Int_t i=0;i<160;++i) {
448 AliTPCclusterMI *c=seed->GetClusterPointer(i);
449 vecX[i]=0;
450 vecY[i]=0;
451 vecZ[i]=0;
452 if (!c) continue;
453 AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
454 if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
455 vecClY[i] = c->GetY();
456 vecClZ[i] = c->GetZ();
457 cl=*c;
458 const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
459 if (c->GetDetector()==s1) ++count1;
460 if (c->GetDetector()==s2) ++count2;
461 Double_t gxyz[3],xyz[3];
462 t1.GetXYZ(gxyz);
463 Float_t bz = AliTracker::GetBz(gxyz);
464 par->GetYAt(c->GetX(), bz, xyz[1]);
465 par->GetZAt(c->GetX(), bz, xyz[2]);
466 vecX[i] = c->GetX();
467 vecY[i]= xyz[1];
468 vecZ[i]= xyz[2];
469 }
470 //
471 //
472 if (fStreamLevel>5){
473 //
474 // huge output - cluster residuals to be investigated
475 //
476 TTreeSRedirector *cstream = GetDebugStreamer();
6a77c9f1 477 //AliTPCseed * t = (AliTPCseed*) seed;
967eae0d 478 //AliExternalTrackParam *p0 = &((AliExternalTrackParam&)seed);
479 AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
480 AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
1c1a1176 481 /*
482
483 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");
484
485 */
486
967eae0d 487 if (cstream){
488 (*cstream)<<"Track"<<
108953e9 489 "run="<<fRun<< // run number
490 "event="<<fEvent<< // event number
491 "time="<<fTime<< // time stamp of event
492 "trigger="<<fTrigger<< // trigger
493 "mag="<<fMagF<< // magnetic field
967eae0d 494 "Cl.="<<&arrCl<<
495 //"tp0.="<<p0<<
496 "tp1.="<<p1<<
497 "tp2.="<<p2<<
498 "vtX.="<<&vecX<<
499 "vtY.="<<&vecY<<
500 "vtZ.="<<&vecZ<<
501 "vcY.="<<&vecClY<<
502 "vcZ.="<<&vecClZ<<
503 "s1="<<s1<<
504 "s2="<<s2<<
505 "c1="<<count1<<
506 "c2="<<count2<<
507 "\n";
508 }
509 }
510}
511
512
513
514
7eaa723e 515void AliTPCcalibAlign::Process12(const Double_t *t1,
516 const Double_t *t2,
9318a5b4 517 TLinearFitter *fitter) {
518 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
519 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
520 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
521 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
522 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
523 //
f8a2dcfb 524 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
525 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
526 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
527
528
529
7eaa723e 530 const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
531 const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 532
533 // TODO:
7eaa723e 534 Double_t sy = 0.1;
535 Double_t sz = 0.1;
536 Double_t sdydx = 0.001;
537 Double_t sdzdx = 0.001;
9318a5b4 538
539 Double_t p[12];
540 Double_t value;
541
542 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
543 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
544 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
545 for (Int_t i=0; i<12;i++) p[i]=0.;
546 p[3+0] = x1; // a10
547 p[3+1] = y1; // a11
548 p[3+2] = z1; // a12
549 p[9+1] = 1.; // a13
550 p[0+1] = y1*dydx2; // a01
551 p[0+2] = z1*dydx2; // a02
552 p[9+0] = dydx2; // a03
553 value = y2;
554 fitter->AddPoint(p,value,sy);
555
556 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
557 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
558 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
559 for (Int_t i=0; i<12;i++) p[i]=0.;
560 p[6+0] = x1; // a20
561 p[6+1] = y1; // a21
562 p[6+2] = z1; // a22
563 p[9+2] = 1.; // a23
564 p[0+1] = y1*dzdx2; // a01
565 p[0+2] = z1*dzdx2; // a02
566 p[9+0] = dzdx2; // a03
567 value = z2;
568 fitter->AddPoint(p,value,sz);
569
570 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
571 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
572 for (Int_t i=0; i<12;i++) p[i]=0.;
573 p[3+0] = 1.; // a10
574 p[3+1] = dydx1; // a11
575 p[3+2] = dzdx1; // a12
576 p[0+0] = -dydx2; // a00
577 p[0+1] = -dydx1*dydx2; // a01
578 p[0+2] = -dzdx1*dydx2; // a02
579 value = 0.;
580 fitter->AddPoint(p,value,sdydx);
581
582 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
583 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
584 for (Int_t i=0; i<12;i++) p[i]=0.;
585 p[6+0] = 1; // a20
586 p[6+1] = dydx1; // a21
587 p[6+2] = dzdx1; // a22
588 p[0+0] = -dzdx2; // a00
589 p[0+1] = -dydx1*dzdx2; // a01
590 p[0+2] = -dzdx1*dzdx2; // a02
591 value = 0.;
592 fitter->AddPoint(p,value,sdzdx);
593}
594
595void AliTPCcalibAlign::Process9(Double_t *t1,
596 Double_t *t2,
597 TLinearFitter *fitter) {
598 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
599 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
600 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
601 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
602 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
603 //
f8a2dcfb 604 // a00 a01 a02 a03 1 p[0] p[1] p[6]
605 // a10 a11 a12 a13 ==> p[2] 1 p[3] p[7]
606 // a20 a21 a21 a23 p[4] p[5] 1 p[8]
607
608
9318a5b4 609 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 610 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 611
612 // TODO:
f8a2dcfb 613 Double_t sy = 0.1;
614 Double_t sz = 0.1;
615 Double_t sdydx = 0.001;
616 Double_t sdzdx = 0.001;
617 //
9318a5b4 618 Double_t p[12];
619 Double_t value;
620
621 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
622 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
f8a2dcfb 623 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
9318a5b4 624 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 625 p[2] += x1; // a10
626 //p[] +=1; // a11
627 p[3] += z1; // a12
628 p[7] += 1; // a13
629 p[0] += y1*dydx2; // a01
630 p[1] += z1*dydx2; // a02
631 p[6] += dydx2; // a03
632 value = y2-y1; //-a11
9318a5b4 633 fitter->AddPoint(p,value,sy);
f8a2dcfb 634 //
9318a5b4 635 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
636 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
f8a2dcfb 637 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
9318a5b4 638 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 639 p[4] += x1; // a20
640 p[5] += y1; // a21
641 //p[] += z1; // a22
642 p[8] += 1.; // a23
643 p[0] += y1*dzdx2; // a01
644 p[1] += z1*dzdx2; // a02
645 p[6] += dzdx2; // a03
646 value = z2-z1; //-a22
9318a5b4 647 fitter->AddPoint(p,value,sz);
648
649 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
650 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
651 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 652 p[2] += 1.; // a10
653 //p[] += dydx1; // a11
654 p[3] += dzdx1; // a12
655 //p[] += -dydx2; // a00
656 p[0] += -dydx1*dydx2; // a01
657 p[1] += -dzdx1*dydx2; // a02
658 value = -dydx1+dydx2; // -a11 + a00
9318a5b4 659 fitter->AddPoint(p,value,sdydx);
660
661 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
662 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
663 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 664 p[4] += 1; // a20
665 p[5] += dydx1; // a21
666 //p[] += dzdx1; // a22
667 //p[] += -dzdx2; // a00
668 p[0] += -dydx1*dzdx2; // a01
669 p[1] += -dzdx1*dzdx2; // a02
670 value = -dzdx1+dzdx2; // -a22 + a00
9318a5b4 671 fitter->AddPoint(p,value,sdzdx);
672}
673
674void AliTPCcalibAlign::Process6(Double_t *t1,
675 Double_t *t2,
676 TLinearFitter *fitter) {
677 // x2 = 1 *x1 +-a01*y1 + 0 +a03
678 // y2 = a01*x1 + 1 *y1 + 0 +a13
679 // z2 = a20*x1 + a21*y1 + 1 *z1 +a23
680 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
681 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
682 //
f8a2dcfb 683 // a00 a01 a02 a03 1 -p[0] 0 p[3]
684 // a10 a11 a12 a13 ==> p[0] 1 0 p[4]
685 // a20 a21 a21 a23 p[1] p[2] 1 p[5]
686
9318a5b4 687 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 688 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 689
690 // TODO:
7eaa723e 691 Double_t sy = 0.1;
692 Double_t sz = 0.1;
693 Double_t sdydx = 0.001;
694 Double_t sdzdx = 0.001;
9318a5b4 695
696 Double_t p[12];
697 Double_t value;
f8a2dcfb 698 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
699 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
9318a5b4 700 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
701 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 702 p[0] += x1; // a10
703 //p[] +=1; // a11
704 //p[] += z1; // a12
705 p[4] += 1; // a13
706 p[0] += -y1*dydx2; // a01
707 //p[] += z1*dydx2; // a02
708 p[3] += dydx2; // a03
709 value = y2-y1; //-a11
9318a5b4 710 fitter->AddPoint(p,value,sy);
f8a2dcfb 711 //
712 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
713 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
714 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
9318a5b4 715 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 716 p[1] += x1; // a20
717 p[2] += y1; // a21
718 //p[] += z1; // a22
719 p[5] += 1.; // a23
720 p[0] += -y1*dzdx2; // a01
721 //p[] += z1*dzdx2; // a02
722 p[3] += dzdx2; // a03
723 value = z2-z1; //-a22
9318a5b4 724 fitter->AddPoint(p,value,sz);
725
726 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
727 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
728 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 729 p[0] += 1.; // a10
730 //p[] += dydx1; // a11
731 //p[] += dzdx1; // a12
732 //p[] += -dydx2; // a00
733 p[0] += dydx1*dydx2; // a01
734 //p[] += -dzdx1*dydx2; // a02
735 value = -dydx1+dydx2; // -a11 + a00
9318a5b4 736 fitter->AddPoint(p,value,sdydx);
737
738 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
739 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
740 for (Int_t i=0; i<12;i++) p[i]=0.;
f8a2dcfb 741 p[1] += 1; // a20
742 p[2] += dydx1; // a21
743 //p[] += dzdx1; // a22
744 //p[] += -dzdx2; // a00
745 p[0] += dydx1*dzdx2; // a01
746 //p[] += -dzdx1*dzdx2; // a02
747 value = -dzdx1+dzdx2; // -a22 + a00
9318a5b4 748 fitter->AddPoint(p,value,sdzdx);
749}
750
7eaa723e 751
752
753
754void AliTPCcalibAlign::EvalFitters() {
755 //
756 // Analyze function
757 //
758 // Perform the fitting using linear fitters
759 //
760 Int_t kMinPoints =50;
9318a5b4 761 TLinearFitter *f;
7eaa723e 762 TFile fff("alignDebug.root","recreate");
9318a5b4 763 for (Int_t s1=0;s1<72;++s1)
7eaa723e 764 for (Int_t s2=0;s2<72;++s2){
e81dc112 765 if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
766 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
0ebabeb6 767 if (f->Eval()!=0) {
9318a5b4 768 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
7eaa723e 769 f->Write(Form("f12_%d_%d",s1,s2));
770 }else{
771 f->Write(Form("f12_%d_%d",s1,s2));
9318a5b4 772 }
773 }
e81dc112 774 if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
775 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
0ebabeb6 776 if (f->Eval()!=0) {
7eaa723e 777 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
778 }else{
779 f->Write(Form("f9_%d_%d",s1,s2));
780 }
781 }
e81dc112 782 if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
783 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
972cf6f2 784 if (f->Eval()!=0) {
7eaa723e 785 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
786 }else{
787 f->Write(Form("f6_%d_%d",s1,s2));
788 }
789 }
790 }
0ebabeb6 791 this->Write("align");
9318a5b4 792 /*
793
794 fitter->Eval();
795 fitter->Eval();
796 chi212 = align->GetChisquare()/(4.*entries);
797
798 TMatrixD mat(13,13);
799 TVectorD par(13);
800 align->GetParameters(par);
801 align->GetCovarianceMatrix(mat);
802
803 //
804 //
805 for (Int_t i=0; i<12;i++){
806 palign12(i)= par(i+1);
807 for (Int_t j=0; j<12;j++){
808 pcovar12(i,j) = mat(i+1,j+1);
809 pcovar12(i,j) *= chi212;
810 }
811 }
812 //
813 for (Int_t i=0; i<12;i++){
814 psigma12(i) = TMath::Sqrt(pcovar12(i,i));
815 palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i));
816 for (Int_t j=0; j<12;j++){
817 pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j));
818 }
819 }
820 */
821}
822
972cf6f2 823TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
824 //
825 // get or make fitter - general linear transformation
826 //
e81dc112 827 static Int_t counter12=0;
828 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 829 TLinearFitter * fitter = GetFitter12(s1,s2);
830 if (fitter) return fitter;
e81dc112 831 // 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]");
832 fitter =new TLinearFitter(&f12,"");
bb6bc8f6 833 fitter->StoreData(kTRUE);
e81dc112 834 fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
835 counter12++;
836 if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter12<<endl;
972cf6f2 837 return fitter;
838}
839
840TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
841 //
842 //get or make fitter - general linear transformation - no scaling
843 //
e81dc112 844 static Int_t counter9=0;
845 static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
972cf6f2 846 TLinearFitter * fitter = GetFitter9(s1,s2);
847 if (fitter) return fitter;
e81dc112 848 // fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
849 fitter =new TLinearFitter(&f9,"");
bb6bc8f6 850 fitter->StoreData(kTRUE);
972cf6f2 851 fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
e81dc112 852 counter9++;
853 if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter9<<endl;
972cf6f2 854 return fitter;
855}
856
857TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
858 //
859 // get or make fitter - 6 paramater linear tranformation
860 // - no scaling
861 // - rotation x-y
862 // - tilting x-z, y-z
e81dc112 863 static Int_t counter6=0;
864 static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
972cf6f2 865 TLinearFitter * fitter = GetFitter6(s1,s2);
866 if (fitter) return fitter;
e81dc112 867 // fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
868 fitter=new TLinearFitter(&f6,"");
bb6bc8f6 869 fitter->StoreData(kTRUE);
972cf6f2 870 fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
e81dc112 871 counter6++;
872 if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<" : "<<counter6<<endl;
972cf6f2 873 return fitter;
874}
875
876
877
878
879
9318a5b4 880Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
972cf6f2 881 //
882 // GetTransformation matrix - 12 paramaters - generael linear transformation
883 //
9318a5b4 884 if (!GetFitter12(s1,s2))
885 return false;
886 else {
887 TVectorD p(12);
9318a5b4 888 GetFitter12(s1,s2)->GetParameters(p);
9318a5b4 889 a.ResizeTo(4,4);
972cf6f2 890 a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
891 a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
892 a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
893 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
9318a5b4 894 return true;
895 }
896}
897
898Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
972cf6f2 899 //
900 // GetTransformation matrix - 9 paramaters - general linear transformation
901 // No scaling
902 //
9318a5b4 903 if (!GetFitter9(s1,s2))
904 return false;
905 else {
906 TVectorD p(9);
907 GetFitter9(s1,s2)->GetParameters(p);
908 a.ResizeTo(4,4);
f8a2dcfb 909 a[0][0]=1; a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
910 a[1][0]=p[2]; a[1][1]=1; a[1][2]=p[3]; a[1][3]=p[7];
911 a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1; a[2][3]=p[8];
972cf6f2 912 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
9318a5b4 913 return true;
914 }
915}
916
917Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
972cf6f2 918 //
919 // GetTransformation matrix - 6 paramaters
920 // 3 translation
921 // 1 rotation -x-y
922 // 2 tilting x-z y-z
9318a5b4 923 if (!GetFitter6(s1,s2))
924 return false;
925 else {
926 TVectorD p(6);
9318a5b4 927 GetFitter6(s1,s2)->GetParameters(p);
9318a5b4 928 a.ResizeTo(4,4);
f8a2dcfb 929 a[0][0]=1; a[0][1]=-p[0];a[0][2]=0; a[0][3]=p[3];
930 a[1][0]=p[0]; a[1][1]=1; a[1][2]=0; a[1][3]=p[4];
931 a[2][0]=p[1]; a[2][1]=p[2]; a[2][2]=1; a[2][3]=p[5];
972cf6f2 932 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
9318a5b4 933 return true;
934 }
935}
972cf6f2 936
937void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
938 const AliExternalTrackParam &tp2,
939 Int_t s1,Int_t s2) {
940 //
941 // Fill residual histograms
8b3c60d8 942 // Innner-Outer
943 // Left right - x-y
944 // A-C side
945 if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0) {
972cf6f2 946 GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
947 GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
948 GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
949 GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
bb6bc8f6 950 //
951 GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
952 GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
953 GetHisto(kYz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetY()-tp2.GetY());
954 GetHisto(kZz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetZ()-tp2.GetZ());
955 //
956 GetHisto(kYPhi,s1,s2,kTRUE)->Fill(tp1.GetSnp(),tp1.GetY()-tp2.GetY());
957 GetHisto(kZTheta,s1,s2,kTRUE)->Fill(tp1.GetTgl(),tp1.GetZ()-tp2.GetZ());
958
959
972cf6f2 960 }
961}
962
963
964
965TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
966{
967 //
968 // return specified residual histogram - it is only QA
969 // if force specified the histogram and given histogram is not existing
970 // new histogram is created
971 //
972 if (GetIndex(s1,s2)>=72*72) return 0;
973 TObjArray *histoArray=0;
974 switch (type) {
975 case kY:
976 histoArray = &fDyHistArray; break;
977 case kZ:
978 histoArray = &fDzHistArray; break;
979 case kPhi:
980 histoArray = &fDphiHistArray; break;
981 case kTheta:
982 histoArray = &fDthetaHistArray; break;
bb6bc8f6 983 case kYPhi:
984 histoArray = &fDyPhiHistArray; break;
985 case kZTheta:
986 histoArray = &fDzThetaHistArray; break;
987 case kYz:
988 histoArray = &fDyZHistArray; break;
989 case kZz:
990 histoArray = &fDzZHistArray; break;
991 case kPhiZ:
992 histoArray = &fDphiZHistArray; break;
993 case kThetaZ:
994 histoArray = &fDthetaZHistArray; break;
972cf6f2 995 }
996 TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
997 if (histo) return histo;
998 if (force==kFALSE) return 0;
999 //
1000 stringstream name;
1001 stringstream title;
1002 switch (type) {
1003 case kY:
1004 name<<"hist_y_"<<s1<<"_"<<s2;
1005 title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1006 histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1007 break;
1008 case kZ:
1009 name<<"hist_z_"<<s1<<"_"<<s2;
1010 title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1011 histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1012 break;
1013 case kPhi:
1014 name<<"hist_phi_"<<s1<<"_"<<s2;
1015 title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1016 histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1017 break;
1018 case kTheta:
1019 name<<"hist_theta_"<<s1<<"_"<<s2;
1020 title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1021 histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1022 break;
bb6bc8f6 1023 //
1024 //
1025 case kYPhi:
1026 name<<"hist_yphi_"<<s1<<"_"<<s2;
1027 title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1028 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,128,-0.3,0.3); // +/- 3 mm
1029 break;
1030 case kZTheta:
1031 name<<"hist_ztheta_"<<s1<<"_"<<s2;
1032 title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1033 histo = new TH2F(name.str().c_str(),title.str().c_str(),128,20,-1,1,-0.3,0.3); // +/- 3 mm
1034 break;
1035 //
1036 //
1037 //
1038 case kYz:
1039 name<<"hist_yz_"<<s1<<"_"<<s2;
1040 title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1041 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1042 break;
1043 case kZz:
1044 name<<"hist_zz_"<<s1<<"_"<<s2;
1045 title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1046 histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1047 break;
1048 case kPhiZ:
1049 name<<"hist_phiz_"<<s1<<"_"<<s2;
1050 title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1051 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1052 break;
1053 case kThetaZ:
1054 name<<"hist_thetaz_"<<s1<<"_"<<s2;
1055 title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1056 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1057 break;
1058
1059
972cf6f2 1060 }
1061 histo->SetDirectory(0);
1062 histoArray->AddAt(histo,GetIndex(s1,s2));
1063 return histo;
1064}
8b3c60d8 1065
1066TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
1067 Int_t i0, Int_t i1, FitType type)
1068{
1069 //
1070 //
1071 //
1072 TMatrixD mat;
6a77c9f1 1073 //TObjArray *fitArray=0;
8b3c60d8 1074 Double_t xsec[1000];
1075 Double_t ysec[1000];
1076 Int_t npoints=0;
1077 for (Int_t isec = sec0; isec<=sec1; isec++){
1078 Int_t isec2 = (isec+dsec)%72;
1079 switch (type) {
1080 case k6:
1081 GetTransformation6(isec,isec2,mat);break;
1082 case k9:
1083 GetTransformation9(isec,isec2,mat);break;
1084 case k12:
1085 GetTransformation12(isec,isec2,mat);break;
1086 }
1087 xsec[npoints]=isec;
1088 ysec[npoints]=mat(i0,i1);
1089 ++npoints;
1090 }
1091 TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1092 Char_t name[1000];
1093 sprintf(name,"Mat[%d,%d] Type=%d",i0,i1,type);
1094 gr->SetName(name);
1095 return gr;
1096}
1097
1098void AliTPCcalibAlign::MakeTree(const char *fname){
1099 //
1100 // make tree with alignment cosntant -
1101 // For QA visualization
1102 //
ae0ac7be 1103 /*
1104 TFile f("CalibObjects.root");
1105 TObjArray *array = (TObjArray*)f.Get("TPCCalib");
1106 AliTPCcalibAlign *alignTPC = (AliTPCcalibAlign *)array->At(0);
1107 alignTPC->MakeTree("alignTree.root");
1108 TFile falign("alignTree.root");
1109 Align->Draw("dy")
1110 */
8b3c60d8 1111 const Int_t kMinPoints=50;
1112 TTreeSRedirector cstream(fname);
1113 for (Int_t s1=0;s1<72;++s1)
1114 for (Int_t s2=0;s2<72;++s2){
1115 if (fPoints[GetIndex(s1,s2)]<kMinPoints) continue;
1116 TMatrixD m6;
1117 TMatrixD m9;
1118 TMatrixD m12;
1119 GetTransformation6(s1,s2,m6);
1120 GetTransformation9(s1,s2,m9);
1121 GetTransformation12(s1,s2,m12);
1122 Double_t dy=0, dz=0, dphi=0,dtheta=0;
1123 Double_t sy=0, sz=0, sphi=0,stheta=0;
1124 Double_t ny=0, nz=0, nphi=0,ntheta=0;
1125 TH1 * his=0;
1126 his = GetHisto(kY,s1,s2);
1127 if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1128 his = GetHisto(kZ,s1,s2);
1129 if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1130 his = GetHisto(kPhi,s1,s2);
1131 if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
1132 his = GetHisto(kTheta,s1,s2);
1133 if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
1134 //
1135 cstream<<"Align"<<
1136 "s1="<<s1<< // reference sector
1137 "s2="<<s2<< // sector to align
1138 "m6.="<<&m6<< // tranformation matrix
1139 "m9.="<<&m9<< //
1140 "m12.="<<&m12<<
967eae0d 1141 // histograms mean RMS and entries
1142 "dy="<<dy<<
8b3c60d8 1143 "sy="<<sy<<
1144 "ny="<<ny<<
1145 "dz="<<dz<<
1146 "sz="<<sz<<
1147 "nz="<<nz<<
1148 "dphi="<<dphi<<
1149 "sphi="<<sphi<<
1150 "nphi="<<nphi<<
1151 "dtheta="<<dtheta<<
1152 "stheta="<<stheta<<
1153 "ntheta="<<ntheta<<
1154 "\n";
1155 }
1156
1157}
ae0ac7be 1158
1159
1160//_____________________________________________________________________
1161Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
1162 //
1163 // merge function
1164 //
1165 if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
1166 if (!list)
1167 return 0;
1168 if (list->IsEmpty())
1169 return 1;
1170
1171 TIterator* iter = list->MakeIterator();
1172 TObject* obj = 0;
1173 iter->Reset();
1174 Int_t count=0;
1175 while((obj = iter->Next()) != 0)
1176 {
1177 AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
1178 if (entry == 0) continue;
1179 Add(entry);
1180 count++;
1181 }
1182 return count;
1183}
1184
1185
1186void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
1187 //
bb6bc8f6 1188 // Add entry - used for merging of compoents
ae0ac7be 1189 //
ae0ac7be 1190 for (Int_t i=0; i<72;i++){
1191 for (Int_t j=0; j<72;j++){
bb6bc8f6 1192 if (align->fPoints[GetIndex(i,j)]==0) continue;
ae0ac7be 1193 fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
ae0ac7be 1194 //
ae0ac7be 1195 //
ae0ac7be 1196 //
bb6bc8f6 1197 for (Int_t itype=0; itype<10; itype++){
1198 TH1 * his0=0, *his1=0;
1199 his0 = GetHisto((HistoType)itype,i,j);
1200 his1 = align->GetHisto((HistoType)itype,i,j);
1201 if (his1){
1202 if (his0) his0->Add(his1);
1203 else {
1204 his0 = GetHisto(kY,i,j,kTRUE);
1205 his0->Add(his1);
1206 }
1207 }
1208 }
1209 // //
1210// // dy
1211// TH1* hdy0 = GetHisto(kY,i,j);
1212// TH1* hdy1 = align->GetHisto(kY,i,j);
1213// if (hdy1){
1214// if (hdy0) hdy0->Add(hdy1);
1215// else {
1216// hdy0 = GetHisto(kY,i,j,kTRUE);
1217// hdy0->Add(hdy1);
1218// }
1219// }
1220// //
1221// // dz
1222// TH1* hdz0 = GetHisto(kZ,i,j);
1223// TH1* hdz1 = align->GetHisto(kZ,i,j);
1224// if (hdz1){
1225// if (hdz0) hdz0->Add(hdz1);
1226// else {
1227// hdz0 = GetHisto(kZ,i,j,kTRUE);
1228// hdz0->Add(hdz1);
1229// }
1230// }
1231// //
1232// // dphi
1233// TH1* hdphi0 = GetHisto(kPhi,i,j);
1234// TH1* hdphi1 = align->GetHisto(kPhi,i,j);
1235// if (hdphi1){
1236// if (hdphi0) hdphi0->Add(hdphi1);
1237// else {
1238// hdphi0 = GetHisto(kPhi,i,j,kTRUE);
1239// hdphi0->Add(hdphi1);
1240// }
1241// }
1242// //
1243// // dtheta
1244// TH1* hdTheta0 = GetHisto(kTheta,i,j);
1245// TH1* hdTheta1 = align->GetHisto(kTheta,i,j);
1246// if (hdTheta1){
1247// if (hdTheta0) hdTheta0->Add(hdTheta1);
1248// else {
1249// hdTheta0 = GetHisto(kTheta,i,j,kTRUE);
1250// hdTheta0->Add(hdTheta1);
1251// }
1252// }
ae0ac7be 1253 }
1254 }
1255 TLinearFitter *f0=0;
1256 TLinearFitter *f1=0;
1257 for (Int_t i=0; i<72;i++){
1258 for (Int_t j=0; j<72;j++){
bb6bc8f6 1259 if (align->fPoints[GetIndex(i,j)]==0) continue;
ae0ac7be 1260 //
1261 // fitter12
1262 f0 = GetFitter12(i,j);
bb6bc8f6 1263 f1 = align->GetFitter12(i,j);
ae0ac7be 1264 if (f1){
1265 if (f0) f0->Add(f1);
1266 else {
1267 f0 = GetOrMakeFitter12(i,j);
1268 f0->Add(f1);
1269 }
1270 }
1271 //
1272 // fitter9
1273 f0 = GetFitter9(i,j);
bb6bc8f6 1274 f1 = align->GetFitter9(i,j);
ae0ac7be 1275 if (f1){
1276 if (f0) f0->Add(f1);
1277 else {
1278 f0 = GetOrMakeFitter9(i,j);
1279 f0->Add(f1);
1280 }
1281 }
1282 f0 = GetFitter6(i,j);
bb6bc8f6 1283 f1 = align->GetFitter6(i,j);
ae0ac7be 1284 if (f1){
1285 if (f0) f0->Add(f1);
1286 else {
1287 f0 = GetOrMakeFitter6(i,j);
bb6bc8f6 1288 f0->Add(f1);
ae0ac7be 1289 }
1290 }
1291 }
1292 }
1293}
108953e9 1294
1295
1296
1297
1298
1299/*
1300
1301
1302gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1303gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1304AliXRDPROOFtoolkit tool;
1305TChain * chainTr = tool.MakeChain("align.txt","Track",0,10200);
1306chainTr->Lookup();
1307
bb6bc8f6 1308TCut c1pt("abs((tp1.fP[4]+tp2.fP[4])*0.5)<3"); // pt cut - OK
1309TCut cdy("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1310TCut cdz("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1311TCut cdphi("abs(tp1.fP[2]-tp2.fP[2])<0.02");
1312TCut cdt("abs(tp1.fP[3]-tp2.fP[3])<0.02");
1313TCut cd1pt("abs(tp1.fP[4]-tp2.fP[4])<0.3"); // delta 1/pt cut -OK
1314//
1315//
1316TCut cutA = c1pt+cdy+cdz+cdphi+cdt+cd1pt;
1317chainTr->Draw(">>listEL",acut,"entryList");
1318TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1319chainTr->SetEntryList(elist);
1320
1321
1322
1323
108953e9 1324
1325
1326TCut cutS("s1%36==s2%36");
1327
1328TCut cutN("c1>32&&c2>60");
1329TCut cutC0("sqrt(tp2.fC[0])<1");
1330
1331TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.4");
1332TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.01");
1333TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.01");
1334TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
1335TCut cutP=cutP0+cutP2+cutP3+cutP4+cutC0;
1336
1337TCut cutX("abs(tp2.fX-133.6)<2");
1338
1339TCut cutA = cutP+cutN;
1340
1341
1342TCut cutY("abs(vcY.fElements-vtY.fElements)<0.3&&vcY.fElements!=0")
1343TCut cutZ("abs(vcZ.fElements-vtZ.fElements)<0.3&&vcZ.fElements!=0")
1344
1345*/