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