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