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