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