1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // Class to make a internal alignemnt of TPC chambers //
21 // Requierements - Warnings:
22 // 1. Before using this componenent the magnetic filed has to be set properly //
23 // 2. The systematic effects - unlinearities has to be understood
25 // If systematic and unlinearities are not under control
26 // the alignment is just effective alignment. Not second order corrction
29 // The histograming of the edge effects and unlineratities integral part
30 // of the component (currently only in debug stream)
32 // 3 general type of linear transformation investigated (see bellow)
34 // By default only 6 parameter alignment to be used - other just for QA purposes
36 // Different linear tranformation investigated
37 // 12 parameters - arbitrary linear transformation
38 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
39 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
40 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
42 // 9 parameters - scaling fixed to 1
43 // a00 a01 a02 a03 1 p[0] p[1] p[6]
44 // a10 a11 a12 a13 ==> p[2] 1 p[3] p[7]
45 // a20 a21 a22 a23 p[4] p[5] 1 p[8]
47 // 6 parameters - x-y rotation x-z, y-z tiliting
48 // a00 a01 a02 a03 1 -p[0] 0 p[3]
49 // a10 a11 a12 a13 ==> p[0] 1 0 p[4]
50 // a20 a21 a22 a23 p[1] p[2] 1 p[5]
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
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
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);
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");
80 .x $ALICE_ROOT/macros/loadlibsREC.C
82 gSystem->Load("$ROOTSYS/lib/libXrdClient.so");
83 gSystem->Load("libProof");
84 gSystem->Load("libANALYSIS");
85 gSystem->Load("libSTAT");
86 gSystem->Load("libTPCcalib");
89 TFile fcalib("CalibObjects.root");
90 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
92 AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
96 align->MakeTree("alignTree.root");
97 TFile falignTree("alignTree.root");
98 TTree * treeAlign = (TTree*)falignTree.Get("Align");
106 #include "TLinearFitter.h"
107 #include "AliTPCcalibAlign.h"
108 #include "AliTPCROC.h"
109 #include "AliTPCPointCorrection.h"
110 #include "AliTrackPointArray.h"
112 #include "AliExternalTrackParam.h"
113 #include "AliESDEvent.h"
114 #include "AliESDfriend.h"
115 #include "AliESDtrack.h"
117 #include "AliTPCTracklet.h"
120 #include "THnSparse.h"
122 #include "TVectorD.h"
123 #include "TTreeStream.h"
127 #include "TGraphErrors.h"
128 #include "AliTPCclusterMI.h"
129 #include "AliTPCseed.h"
130 #include "AliTracker.h"
131 #include "TClonesArray.h"
134 #include "TProfile.h"
136 #include "TDatabasePDG.h"
138 #include "TTreeStream.h"
139 #include "Riostream.h"
144 AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0;
145 Double_t AliTPCcalibAlign::fgkMergeEntriesCut=10000000.; //10**7 tracks
146 ClassImp(AliTPCcalibAlign)
151 AliTPCcalibAlign* AliTPCcalibAlign::Instance()
154 // Singleton implementation
155 // Returns an instance of this class, it is created if neccessary
157 if (fgInstance == 0){
158 fgInstance = new AliTPCcalibAlign();
166 AliTPCcalibAlign::AliTPCcalibAlign()
168 fDphiHistArray(72*72),
169 fDthetaHistArray(72*72),
173 fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi
174 fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta
175 fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz
176 fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz
177 fDyZHistArray(72*72), // array of residual histograms y -kYz
178 fDzZHistArray(72*72), // array of residual histograms z -kZz
179 fFitterArray12(72*72),
180 fFitterArray9(72*72),
181 fFitterArray6(72*72),
182 fMatrixArray12(72*72),
183 fMatrixArray9(72*72),
184 fMatrixArray6(72*72),
185 fCombinedMatrixArray6(72),
190 fArraySectorIntParam(36), // array of sector alignment parameters
191 fArraySectorIntCovar(36), // array of sector alignment covariances
193 // Kalman filter for global alignment
195 fSectorParamA(0), // Kalman parameter for A side
196 fSectorCovarA(0), // Kalman covariance for A side
197 fSectorParamC(0), // Kalman parameter for A side
198 fSectorCovarC(0), // Kalman covariance for A side
199 fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
204 for (Int_t i=0;i<72*72;++i) {
207 AliTPCROC * roc = AliTPCROC::Instance();
208 fXquadrant = roc->GetPadRowRadii(36,53);
209 fXmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
210 fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
211 fClusterDelta[0]=0; // cluster residuals - Y
212 fClusterDelta[1]=0; // cluster residuals - Z
215 fTrackletDelta[0]=0; // tracklet residuals
216 fTrackletDelta[1]=0; // tracklet residuals
217 fTrackletDelta[2]=0; // tracklet residuals
218 fTrackletDelta[3]=0; // tracklet residuals
221 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
223 fDphiHistArray(72*72),
224 fDthetaHistArray(72*72),
227 fDyPhiHistArray(72*72), // array of residual histograms y -kYPhi
228 fDzThetaHistArray(72*72), // array of residual histograms z-z -kZTheta
229 fDphiZHistArray(72*72), // array of residual histograms phi -kPhiz
230 fDthetaZHistArray(72*72), // array of residual histograms theta -kThetaz
231 fDyZHistArray(72*72), // array of residual histograms y -kYz
232 fDzZHistArray(72*72), // array of residual histograms z -kZz //
233 fFitterArray12(72*72),
234 fFitterArray9(72*72),
235 fFitterArray6(72*72),
236 fMatrixArray12(72*72),
237 fMatrixArray9(72*72),
238 fMatrixArray6(72*72),
239 fCombinedMatrixArray6(72),
244 fArraySectorIntParam(36), // array of sector alignment parameters
245 fArraySectorIntCovar(36), // array of sector alignment covariances
247 // Kalman filter for global alignment
249 fSectorParamA(0), // Kalman parameter for A side
250 fSectorCovarA(0), // Kalman covariance for A side
251 fSectorParamC(0), // Kalman parameter for A side
252 fSectorCovarC(0), // Kalman covariance for A side
253 fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
261 for (Int_t i=0;i<72*72;++i) {
264 AliTPCROC * roc = AliTPCROC::Instance();
265 fXquadrant = roc->GetPadRowRadii(36,53);
266 fXmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
267 fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
268 fClusterDelta[0]=0; // cluster residuals
269 fClusterDelta[1]=0; // cluster residuals
271 fTrackletDelta[0]=0; // tracklet residuals
272 fTrackletDelta[1]=0; // tracklet residuals
273 fTrackletDelta[2]=0; // tracklet residuals
274 fTrackletDelta[3]=0; // tracklet residuals
278 AliTPCcalibAlign::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
291 fFitterArray12(align.fFitterArray12),
292 fFitterArray9(align.fFitterArray9),
293 fFitterArray6(align.fFitterArray6),
295 fMatrixArray12(align.fMatrixArray12),
296 fMatrixArray9(align.fMatrixArray9),
297 fMatrixArray6(align.fMatrixArray6),
298 fCombinedMatrixArray6(align.fCombinedMatrixArray6),
299 fNoField(align.fNoField),
301 fXmiddle(align.fXmiddle),
302 fXquadrant(align.fXquadrant),
303 fArraySectorIntParam(align.fArraySectorIntParam), // array of sector alignment parameters
304 fArraySectorIntCovar(align.fArraySectorIntCovar), // array of sector alignment covariances
305 fSectorParamA(0), // Kalman parameter for A side
306 fSectorCovarA(0), // Kalman covariance for A side
307 fSectorParamC(0), // Kalman parameter for A side
308 fSectorCovarC(0), // Kalman covariance for A side
309 fUseInnerOuter(kTRUE)// flag- use Inner Outer sector for left righ alignment
313 // copy constructor - copy also the content
317 const TObjArray *arr1=0;
318 for (Int_t index =0; index<72*72; index++){
319 for (Int_t iarray=0;iarray<10; iarray++){
321 arr0 = &fDyHistArray;
322 arr1 = &align.fDyHistArray;
325 arr0 = &fDzHistArray;
326 arr1 = &align.fDzHistArray;
329 arr0 = &fDphiHistArray;
330 arr1 = &align.fDphiHistArray;
333 arr0 = &fDthetaHistArray;
334 arr1 = &align.fDthetaHistArray;
337 arr0 = &fDyZHistArray;
338 arr1 = &align.fDyZHistArray;
341 arr0 = &fDzZHistArray;
342 arr1 = &align.fDzZHistArray;
345 arr0 = &fDphiZHistArray;
346 arr1 = &align.fDphiZHistArray;
348 if (iarray==kThetaZ){
349 arr0 = &fDthetaZHistArray;
350 arr1 = &align.fDthetaZHistArray;
354 arr0 = &fDyPhiHistArray;
355 arr1 = &align.fDyPhiHistArray;
357 if (iarray==kZTheta){
358 arr0 = &fDzThetaHistArray;
359 arr1 = &align.fDzThetaHistArray;
362 if (arr1->At(index)) {
363 his = (TH1*)arr1->At(index)->Clone();
364 his->SetDirectory(0);
365 arr0->AddAt(his,index);
372 if (align.fSectorParamA){
373 fSectorParamA = (TMatrixD*)align.fSectorParamA->Clone();
374 fSectorParamA = (TMatrixD*)align.fSectorCovarA->Clone();
375 fSectorParamC = (TMatrixD*)align.fSectorParamA->Clone();
376 fSectorParamC = (TMatrixD*)align.fSectorCovarA->Clone();
378 fClusterDelta[0]=0; // cluster residuals
379 fClusterDelta[1]=0; // cluster residuals
381 fTrackletDelta[0]=0; // tracklet residuals
382 fTrackletDelta[1]=0; // tracklet residuals
383 fTrackletDelta[2]=0; // tracklet residuals
384 fTrackletDelta[3]=0; // tracklet residuals
388 AliTPCcalibAlign::~AliTPCcalibAlign() {
392 fDphiHistArray.SetOwner(kTRUE); // array of residual histograms phi -kPhi
393 fDthetaHistArray.SetOwner(kTRUE); // array of residual histograms theta -kTheta
394 fDyHistArray.SetOwner(kTRUE); // array of residual histograms y -kY
395 fDzHistArray.SetOwner(kTRUE); // array of residual histograms z -kZ
397 fDyPhiHistArray.SetOwner(kTRUE); // array of residual histograms y -kYPhi
398 fDzThetaHistArray.SetOwner(kTRUE); // array of residual histograms z-z -kZTheta
400 fDphiZHistArray.SetOwner(kTRUE); // array of residual histograms phi -kPhiz
401 fDthetaZHistArray.SetOwner(kTRUE); // array of residual histograms theta -kThetaz
402 fDyZHistArray.SetOwner(kTRUE); // array of residual histograms y -kYz
403 fDzZHistArray.SetOwner(kTRUE); // array of residual histograms z -kZz
405 fDphiHistArray.Delete(); // array of residual histograms phi -kPhi
406 fDthetaHistArray.Delete(); // array of residual histograms theta -kTheta
407 fDyHistArray.Delete(); // array of residual histograms y -kY
408 fDzHistArray.Delete(); // array of residual histograms z -kZ
410 fDyPhiHistArray.Delete(); // array of residual histograms y -kYPhi
411 fDzThetaHistArray.Delete(); // array of residual histograms z-z -kZTheta
413 fDphiZHistArray.Delete(); // array of residual histograms phi -kPhiz
414 fDthetaZHistArray.Delete(); // array of residual histograms theta -kThetaz
415 fDyZHistArray.Delete(); // array of residual histograms y -kYz
416 fDzZHistArray.Delete(); // array of residual histograms z -kZz
418 fFitterArray12.SetOwner(kTRUE); // array of fitters
419 fFitterArray9.SetOwner(kTRUE); // array of fitters
420 fFitterArray6.SetOwner(kTRUE); // array of fitters
422 fMatrixArray12.SetOwner(kTRUE); // array of transnformtation matrix
423 fMatrixArray9.SetOwner(kTRUE); // array of transnformtation matrix
424 fMatrixArray6.SetOwner(kTRUE); // array of transnformtation matrix
426 fFitterArray12.Delete(); // array of fitters
427 fFitterArray9.Delete(); // array of fitters
428 fFitterArray6.Delete(); // array of fitters
430 fMatrixArray12.Delete(); // array of transnformtation matrix
431 fMatrixArray9.Delete(); // array of transnformtation matrix
432 fMatrixArray6.Delete(); // array of transnformtation matrix
435 fArraySectorIntParam.SetOwner(kTRUE); // array of sector alignment parameters
436 fArraySectorIntCovar.SetOwner(kTRUE); // array of sector alignment covariances
437 fArraySectorIntParam.Delete(); // array of sector alignment parameters
438 fArraySectorIntCovar.Delete(); // array of sector alignment covariances
439 for (Int_t i=0; i<2; i++){
440 delete fClusterDelta[i]; // cluster residuals
443 for (Int_t i=0; i<4; i++){
444 delete fTrackletDelta[i]; // tracklet residuals
450 void AliTPCcalibAlign::Process(AliESDEvent *event) {
452 // Process pairs of cosmic tracks
454 const Double_t kptDownscale=50; // downscale factor for the low pt particels
455 if (!fClusterDelta[0]) MakeResidualHistos();
456 if (!fTrackletDelta[0]) MakeResidualHistosTracklet();
459 ExportTrackPoints(event); // export track points for external calibration
460 const Int_t kMaxTracks =6;
461 const Int_t kminCl = 40;
462 AliESDfriend *eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
463 if (!eESDfriend) return;
464 if (eESDfriend->TestSkipBit()) return;
465 Int_t ntracks=event->GetNumberOfTracks();
472 for (Int_t i0=0;i0<ntracks;++i0) {
473 AliESDtrack *track0 = event->GetTrack(i0);
474 AliESDfriendTrack *friendTrack = 0;
475 TObject *calibObject=0;
476 AliTPCseed *seed0 = 0;
478 friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
479 if (!friendTrack) continue;
480 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
481 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
483 if (!seed0) continue;
484 fCurrentTrack=track0;
485 fCurrentFriendTrack=friendTrack;
488 Double_t scalept= TMath::Min(1./TMath::Abs(track0->GetParameter()[4]),2.);
489 Bool_t isSelected = (TMath::Exp(2*scalept)>kptDownscale*gRandom->Rndm());
490 if (isSelected) ProcessSeed(seed0);
493 // process cosmic pairs
495 if (ntracks>kMaxTracks) return;
497 //select pairs - for alignment
498 for (Int_t i0=0;i0<ntracks;++i0) {
499 AliESDtrack *track0 = event->GetTrack(i0);
500 // if (track0->GetTPCNcls()<kminCl) continue;
501 track0->GetImpactParameters(dca0[0],dca0[1]);
502 // if (TMath::Abs(dca0[0])>30) continue;
504 for (Int_t i1=0;i1<ntracks;++i1) {
505 if (i0==i1) continue;
506 AliESDtrack *track1 = event->GetTrack(i1);
507 // if (track1->GetTPCNcls()<kminCl) continue;
508 track1->GetImpactParameters(dca1[0],dca1[1]);
509 // fast cuts on dca and theta
510 // if (TMath::Abs(dca1[0]+dca0[0])>15) continue;
511 // if (TMath::Abs(dca1[1]-dca0[1])>15) continue;
512 if (TMath::Abs(track0->GetParameter()[3]+track1->GetParameter()[3])>0.1) continue;
514 AliESDfriendTrack *friendTrack = 0;
515 TObject *calibObject=0;
516 AliTPCseed *seed0 = 0,*seed1=0;
518 friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
519 if (!friendTrack) continue;
520 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
521 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
523 friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i1);;
524 if (!friendTrack) continue;
525 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
526 if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
528 if (!seed0) continue;
531 if (!seed1) continue;
532 Int_t nclsectors0[72], nclsectors1[72];
533 for (Int_t isec=0;isec<72;isec++){
537 for (Int_t i=0;i<160;i++){
538 AliTPCclusterMI *c0=seed0->GetClusterPointer(i);
539 AliTPCclusterMI *c1=seed1->GetClusterPointer(i);
540 if (c0) nclsectors0[c0->GetDetector()]+=1;
541 if (c1) nclsectors1[c1->GetDetector()]+=1;
544 for (Int_t isec0=0; isec0<72;isec0++){
545 if (nclsectors0[isec0]<kminCl) continue;
546 for (Int_t isec1=0; isec1<72;isec1++){
547 if (nclsectors1[isec1]<kminCl) continue;
550 Double_t parLine0[10];
551 Double_t parLine1[10];
552 TMatrixD par0(4,1),cov0(4,4),par1(4,1),cov1(4,4);
553 Bool_t useInnerOuter = kFALSE;
554 if (s1%36!=s0%36) useInnerOuter = fUseInnerOuter; // for left - right alignment both sectors refit can be used if specified
555 Int_t nl0 = RefitLinear(seed0,s0, parLine0, s0,par0,cov0,fXIO,useInnerOuter);
556 Int_t nl1 = RefitLinear(seed1,s1, parLine1, s0,par1,cov1,fXIO,useInnerOuter);
557 parLine0[0]=0; // reference frame in IO boundary
559 // if (nl0<kminCl || nl1<kminCl) continue;
563 if (TMath::Min(nl0,nl1)<kminCl) isOK=kFALSE;
564 // apply selection criteria
568 dp0=par0(0,0)-par1(0,0);
569 dp1=par0(1,0)-par1(1,0);
570 dp3=par0(3,0)-par1(3,0);
571 pp0=dp0/TMath::Sqrt(cov0(0,0)+cov1(0,0)+0.1*0.1);
572 pp1=dp1/TMath::Sqrt(cov0(1,1)+cov1(1,1)+0.0015*0.0015);
573 pp3=dp3/TMath::Sqrt(cov0(3,3)+cov1(3,3)+0.0015*0.0015);
575 if (TMath::Abs(dp0)>1.0) isOK=kFALSE;
576 if (TMath::Abs(dp1)>0.02) isOK=kFALSE;
577 if (TMath::Abs(dp3)>0.02) isOK=kFALSE;
578 if (TMath::Abs(pp0)>6) isOK=kFALSE;
579 if (TMath::Abs(pp1)>6) isOK=kFALSE;
580 if (TMath::Abs(pp3)>6) isOK=kFALSE;
583 FillHisto(parLine0,parLine1,s0,s1);
584 ProcessAlign(parLine0,parLine1,s0,s1);
585 UpdateKalman(s0,s1,par0, cov0, par1, cov1);
588 TTreeSRedirector *cstream = GetDebugStreamer();
590 (*cstream)<<"cosmic"<<
609 void AliTPCcalibAlign::ExportTrackPoints(AliESDEvent *event){
611 // Export track points for alignment - calibration
612 // export space points for pairs of tracks if possible
614 AliESDfriend *eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
615 if (!eESDfriend) return;
616 Int_t ntracks=event->GetNumberOfTracks();
617 Int_t kMaxTracks=4; // maximal number of tracks for cosmic pairs
618 Int_t kMinVertexTracks=5; // maximal number of tracks for vertex mesurement
621 const Int_t kminCl = 60;
622 const Int_t kminClSum = 120;
623 //const Double_t kDistY = 5;
624 // const Double_t kDistZ = 40;
625 const Double_t kDistTh = 0.05;
626 const Double_t kDistThS = 0.002;
627 const Double_t kDist1Pt = 0.1;
628 const Double_t kMaxD0 =3; // max distance to the primary vertex
629 const Double_t kMaxD1 =5; // max distance to the primary vertex
630 const AliESDVertex *tpcVertex = 0;
631 // get the primary vertex TPC
632 if (ntracks>kMinVertexTracks) {
633 tpcVertex = event->GetPrimaryVertexSPD();
634 if (tpcVertex->GetNContributors()<kMinVertexTracks) tpcVertex=0;
639 Int_t index0=0,index1=0;
641 for (Int_t i0=0;i0<ntracks;++i0) {
642 AliESDtrack *track0 = event->GetTrack(i0);
643 if (!track0) continue;
644 if ((track0->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
645 if (track0->GetOuterParam()==0) continue;
646 if (track0->GetInnerParam()==0) continue;
647 if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt()-track0->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
648 if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
649 if (TMath::Abs(track0->GetInnerParam()->GetTgl()-track0->GetOuterParam()->GetTgl())>kDistThS) continue;
650 AliESDtrack *track1P = 0;
651 if (track0->GetTPCNcls()<kminCl) continue;
652 track0->GetImpactParameters(dca0[0],dca0[1]);
656 if (ntracks<kMaxTracks) for (Int_t i1=i0+1;i1<ntracks;++i1) {
657 if (i0==i1) continue;
658 AliESDtrack *track1 = event->GetTrack(i1);
659 if (!track1) continue;
660 if ((track1->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
661 if (track1->GetOuterParam()==0) continue;
662 if (track1->GetInnerParam()==0) continue;
663 if (track1->GetTPCNcls()<kminCl) continue;
664 if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt()-track1->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
665 if (TMath::Abs(track1->GetInnerParam()->GetTgl()-track1->GetOuterParam()->GetTgl())>kDistThS) continue;
666 if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
667 //track1->GetImpactParameters(dca1[0],dca1[1]);
668 //if (TMath::Abs(dca1[0]-dca0[0])>kDistY) continue;
669 //if (TMath::Abs(dca1[1]-dca0[1])>kDistZ) continue;
670 if (TMath::Abs(track0->GetTgl()+track1->GetTgl())>kDistTh) continue;
671 if (TMath::Abs(track0->GetSigned1Pt()+track1->GetSigned1Pt())>kDist1Pt) continue;
675 AliESDfriendTrack *friendTrack = 0;
676 TObject *calibObject=0;
677 AliTPCseed *seed0 = 0,*seed1=0;
679 friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index0);;
680 if (!friendTrack) continue;
681 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
682 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
685 friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index1);;
686 if (!friendTrack) continue;
687 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
688 if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
692 Int_t npoints=0, ncont=0;
693 if (seed0) {npoints+=seed0->GetNumberOfClusters(); ncont++;}
694 if (seed1) {npoints+=seed1->GetNumberOfClusters(); ncont++;}
695 if (npoints<kminClSum) continue;
697 AliTrackPointArray array(npoints);
699 Double_t dxyz[3]={0,0,0};
700 Double_t dc[6]={0,0,0};
701 tpcVertex->GetXYZ(dxyz);
702 tpcVertex->GetCovarianceMatrix(dc);
703 Float_t xyz[3]={dxyz[0],dxyz[1],dxyz[2]};
704 Float_t cov[6]={dc[0]+1,dc[1],dc[2]+1,dc[3],dc[4],dc[5]+100.};
705 AliTrackPoint point(xyz,cov,73); // add point to not existing volume
706 Float_t dtpc[2],dcov[3];
707 track0->GetImpactParametersTPC(dtpc,dcov);
708 if (TMath::Abs(dtpc[0])>kMaxD0) continue;
709 if (TMath::Abs(dtpc[1])>kMaxD1) continue;
712 if (seed0) for (Int_t icl = 0; icl<160; icl++){
713 AliTPCclusterMI *cluster=seed0->GetClusterPointer(icl);
714 if (!cluster) continue;
717 cluster->GetGlobalXYZ(xyz);
718 cluster->GetGlobalCov(cov);
719 AliTrackPoint point(xyz,cov,cluster->GetDetector());
720 array.AddPoint(npoints, &point);
721 if (cpoint>=npoints) continue; //shoul not happen
722 array.AddPoint(cpoint, &point);
725 if (seed1) for (Int_t icl = 0; icl<160; icl++){
726 AliTPCclusterMI *cluster=seed1->GetClusterPointer(icl);
727 if (!cluster) continue;
730 cluster->GetGlobalXYZ(xyz);
731 cluster->GetGlobalCov(cov);
732 AliTrackPoint point(xyz,cov,cluster->GetDetector());
733 array.AddPoint(npoints, &point);
734 if (cpoint>=npoints) continue; //shoul not happen
735 array.AddPoint(cpoint, &point);
741 TTreeSRedirector *cstream = GetDebugStreamer();
743 Bool_t isVertex=(tpcVertex)? kTRUE:kFALSE;
744 Double_t tof0=track0->GetTOFsignal();
745 Double_t tof1=(track1P) ? track1P->GetTOFsignal(): 0;
746 static AliExternalTrackParam dummy;
747 AliExternalTrackParam *p0In = &dummy;
748 AliExternalTrackParam *p1In = &dummy;
749 AliExternalTrackParam *p0Out = &dummy;
750 AliExternalTrackParam *p1Out = &dummy;
752 AliESDVertex *pvertex= (tpcVertex)? (AliESDVertex *)tpcVertex: &vdummy;
754 p0In= new AliExternalTrackParam(*track0);
755 p0Out=new AliExternalTrackParam(*(track0->GetOuterParam()));
758 p1In= new AliExternalTrackParam(*track1P);
759 p1Out=new AliExternalTrackParam(*(track1P->GetOuterParam()));
762 (*cstream)<<"trackPoints"<<
763 "run="<<fRun<< // run number
764 "event="<<fEvent<< // event number
765 "time="<<fTime<< // time stamp of event
766 "trigger="<<fTrigger<< // trigger
767 "triggerClass="<<&fTriggerClass<< // trigger
768 "mag="<<fMagF<< // magnetic field
769 "pvertex.="<<pvertex<< // vertex
771 "isVertex="<<isVertex<< // flag is with prim vertex
772 "tof0="<<tof0<< // tof signal 0
773 "tof1="<<tof1<< // tof signal 1
774 "seed0.="<<seed0<< // track info
775 "ntracks="<<ntracks<< // number of tracks
776 "ncont="<<ncont<< // number of contributors
777 "p0In.="<<p0In<< // track parameters0
778 "p1In.="<<p1In<< // track parameters1
779 "p0Out.="<<p0Out<< // track parameters0
780 "p1Out.="<<p1Out<< // track parameters0
790 void AliTPCcalibAlign::ProcessSeed(AliTPCseed *seed) {
794 // make a kalman tracklets out of seed
796 UpdateClusterDeltaField(seed);
798 AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
800 tracklets.SetOwner();
801 Int_t ntracklets = tracklets.GetEntries();
802 if (ntracklets<2) return;
805 for (Int_t i1=0;i1<ntracklets;i1++)
806 for (Int_t i2=0;i2<ntracklets;i2++){
807 if (i1==i2) continue;
808 AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[i1]);
809 AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[i2]);
810 AliExternalTrackParam *common1=0,*common2=0;
811 if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2)){
812 ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
813 UpdateAlignSector(seed,t1->GetSector());
820 void AliTPCcalibAlign::Analyze(){
828 void AliTPCcalibAlign::Terminate(){
830 // Terminate function
831 // call base terminate + Eval of fitters
833 Info("AliTPCcalibAlign","Terminate");
835 AliTPCcalibBase::Terminate();
839 void AliTPCcalibAlign::UpdatePointCorrection(AliTPCPointCorrection * correction){
841 // Update point correction with alignment coefficients
843 for (Int_t isec=0;isec<36;isec++){
844 TMatrixD * matCorr = (TMatrixD*)(correction->fArraySectorIntParam.At(isec));
845 TMatrixD * matAlign = (TMatrixD*)(fArraySectorIntParam.At(isec));
846 TMatrixD * matAlignCovar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
847 if (!matAlign) continue;
849 correction->fArraySectorIntParam.AddAt(matAlign->Clone(),isec);
850 correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec);
853 (*matCorr)+=(*matAlign);
854 correction->fArraySectorIntCovar.AddAt(matAlignCovar->Clone(),isec);
861 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
862 const AliExternalTrackParam &tp2,
863 const AliTPCseed * seed,
866 // Process function to fill fitters
869 Double_t t1[10],t2[10];
870 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
871 Double_t &x2=t2[0], &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
875 Double_t snp1=tp1.GetSnp();
876 dydx1=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
877 Double_t tgl1=tp1.GetTgl();
878 // dz/dx = 1/(cos(theta)*cos(phi))
879 dzdx1=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
883 Double_t snp2=tp2.GetSnp();
884 dydx2=snp2/TMath::Sqrt((1.-snp2)*(1.+snp2));
885 Double_t tgl2=tp2.GetTgl();
886 dzdx2=tgl2/TMath::Sqrt((1.-snp2)*(1.+snp2));
895 t1[6]=TMath::Sqrt(tp1.GetSigmaY2());
896 t1[7]=TMath::Sqrt(tp1.GetSigmaSnp2());
897 t1[8]=TMath::Sqrt(tp1.GetSigmaZ2());
898 t1[9]=TMath::Sqrt(tp1.GetSigmaTgl2());
900 t2[6]=TMath::Sqrt(tp2.GetSigmaY2());
901 t2[7]=TMath::Sqrt(tp2.GetSigmaSnp2());
902 t2[8]=TMath::Sqrt(tp2.GetSigmaZ2());
903 t2[9]=TMath::Sqrt(tp2.GetSigmaTgl2());
907 Double_t parLine1[10];
908 Double_t parLine2[10];
909 TMatrixD par1(4,1),cov1(4,4),par2(4,1),cov2(4,4);
910 Bool_t useInnerOuter = kFALSE;
911 if (s1%36!=s2%36) useInnerOuter = fUseInnerOuter; // for left - right alignment bot sectors refit can be used if specified
912 Int_t nl1 = RefitLinear(seed,s1, parLine1, s1,par1,cov1,tp1.GetX(), useInnerOuter);
913 Int_t nl2 = RefitLinear(seed,s2, parLine2, s1,par2,cov2,tp1.GetX(), useInnerOuter);
914 parLine1[0]=tp1.GetX()-fXIO; // parameters in IROC-OROC boundary
915 parLine2[0]=tp1.GetX()-fXIO; // parameters in IROC-OROC boundary
919 Int_t accept = AcceptTracklet(tp1,tp2);
920 Int_t acceptLinear = AcceptTracklet(parLine1,parLine2);
924 TTreeSRedirector *cstream = GetDebugStreamer();
926 static TVectorD vec1(5);
927 static TVectorD vec2(5);
928 static TVectorD vecL1(9);
929 static TVectorD vecL2(9);
930 vec1.SetElements(t1);
931 vec2.SetElements(t2);
932 vecL1.SetElements(parLine1);
933 vecL2.SetElements(parLine2);
934 AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
935 AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
936 (*cstream)<<"Tracklet"<<
938 "acceptLinear="<<acceptLinear<< // accept linear tracklets
939 "run="<<fRun<< // run number
940 "event="<<fEvent<< // event number
941 "time="<<fTime<< // time stamp of event
942 "trigger="<<fTrigger<< // trigger
943 "triggerClass="<<&fTriggerClass<< // trigger
944 "mag="<<fMagF<< // magnetic field
945 "isOK="<<accept<< // flag - used for alignment
952 "nl1="<<nl1<< // linear fit - n points
953 "nl2="<<nl2<< // linear fit - n points
954 "vl1.="<<&vecL1<< // linear fits
955 "vl2.="<<&vecL2<< // linear fits
959 if (TMath::Abs(fMagF)<0.005){
963 if (nl1>10 && nl2>10 &&(acceptLinear==0)){
964 ProcessDiff(tp1,tp2, seed,s1,s2);
965 if (TMath::Abs(parLine1[2])<0.8 &&TMath::Abs(parLine1[2])<0.8 ){ //angular cut
966 FillHisto(parLine1,parLine2,s1,s2);
967 ProcessAlign(parLine1,parLine2,s1,s2);
968 FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
969 FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
970 //UpdateKalman(s1,s2,par1, cov1, par2, cov2); - OBSOLETE to be removed - 50 % of time here
974 if (accept>0) return;
976 // fill resolution histograms - previous cut included
977 if (TMath::Abs(fMagF)>0.005){
979 // use Kalman if mag field
981 ProcessDiff(tp1,tp2, seed,s1,s2);
982 FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
983 FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
984 FillHisto(t1,t2,s1,s2);
985 ProcessAlign(t1,t2,s1,s2);
989 void AliTPCcalibAlign::ProcessAlign(Double_t * t1,
993 // Do intersector alignment
995 //Process12(t1,t2,GetOrMakeFitter12(s1,s2));
996 //Process9(t1,t2,GetOrMakeFitter9(s1,s2));
997 Process6(t1,t2,GetOrMakeFitter6(s1,s2));
998 ++fPoints[GetIndex(s1,s2)];
1003 Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
1004 const AliExternalTrackParam &p2) const
1008 // Accept pair of tracklets?
1012 TCut cutS0("sqrt(tp2.fC[0]+tp1.fC[0])<0.2");
1013 TCut cutS1("sqrt(tp2.fC[2]+tp1.fC[2])<0.2");
1014 TCut cutS2("sqrt(tp2.fC[5]+tp1.fC[5])<0.01");
1015 TCut cutS3("sqrt(tp2.fC[9]+tp1.fC[9])<0.01");
1016 TCut cutS4("sqrt(tp2.fC[14]+tp1.fC[14])<0.25");
1017 TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
1019 // parameters matching cuts
1020 TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1021 TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1022 TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03");
1023 TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03");
1024 TCut cutP4("abs(tp1.fP[4]-tp2.fP[4])<0.5");
1025 TCut cutPP4("abs(tp1.fP[4]-tp2.fP[4])/sqrt(tp2.fC[14]+tp1.fC[14])<3");
1026 TCut cutP=cutP0+cutP1+cutP2+cutP3+cutP4+cutPP4;
1031 const Double_t *cp1 = p1.GetCovariance();
1032 const Double_t *cp2 = p2.GetCovariance();
1033 if (TMath::Sqrt(cp1[0]+cp2[0])>0.2) reject|=1;;
1034 if (TMath::Sqrt(cp1[2]+cp2[2])>0.2) reject|=2;
1035 if (TMath::Sqrt(cp1[5]+cp2[5])>0.01) reject|=4;
1036 if (TMath::Sqrt(cp1[9]+cp2[9])>0.01) reject|=8;
1037 if (TMath::Sqrt(cp1[14]+cp2[14])>0.2) reject|=16;
1039 //parameters difference
1040 const Double_t *tp1 = p1.GetParameter();
1041 const Double_t *tp2 = p2.GetParameter();
1042 if (TMath::Abs(tp1[0]-tp2[0])>0.6) reject|=32;
1043 if (TMath::Abs(tp1[1]-tp2[1])>0.6) reject|=64;
1044 if (TMath::Abs(tp1[2]-tp2[2])>0.03) reject|=128;
1045 if (TMath::Abs(tp1[3]-tp2[3])>0.03) reject|=526;
1046 if (TMath::Abs(tp1[4]-tp2[4])>0.4) reject|=1024;
1047 if (TMath::Abs(tp1[4]-tp2[4])/TMath::Sqrt(cp1[14]+cp2[14])>4) reject|=2048;
1050 if (TMath::Abs(tp2[1])>235) reject|=2*4096;
1060 Int_t AliTPCcalibAlign::AcceptTracklet(const Double_t *t1, const Double_t *t2) const
1063 // accept tracklet -
1064 // dist cut + 6 sigma cut
1066 Double_t dy = t2[1]-t1[1];
1067 Double_t dphi = t2[2]-t1[2];
1068 Double_t dz = t2[3]-t1[3];
1069 Double_t dtheta = t2[4]-t1[4];
1071 Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]+0.05*0.05);
1072 Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]+0.001*0.001);
1073 Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]+0.05*0.05);
1074 Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]+0.001*0.001);
1077 if (TMath::Abs(dy)>1.) reject|=2;
1078 if (TMath::Abs(dphi)>0.1) reject|=4;
1079 if (TMath::Abs(dz)>1.) reject|=8;
1080 if (TMath::Abs(dtheta)>0.1) reject|=16;
1082 if (TMath::Abs(dy/sy)>6) reject|=32;
1083 if (TMath::Abs(dphi/sdydx)>6) reject|=64;
1084 if (TMath::Abs(dz/sz)>6) reject|=128;
1085 if (TMath::Abs(dtheta/sdzdx)>6) reject|=256;
1090 void AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
1091 const AliExternalTrackParam &t2,
1092 const AliTPCseed *seed,
1096 // Process local residuals function
1101 TVectorD vecClY(160);
1102 TVectorD vecClZ(160);
1103 TClonesArray arrCl("AliTPCclusterMI",160);
1104 arrCl.ExpandCreateFast(160);
1105 Int_t count1=0, count2=0;
1107 for (Int_t i=0;i<160;++i) {
1108 AliTPCclusterMI *c=seed->GetClusterPointer(i);
1113 AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
1114 if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
1115 vecClY[i] = c->GetY();
1116 vecClZ[i] = c->GetZ();
1118 const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
1119 if (c->GetDetector()==s1) ++count1;
1120 if (c->GetDetector()==s2) ++count2;
1121 Double_t gxyz[3],xyz[3];
1123 Float_t bz = AliTracker::GetBz(gxyz);
1124 par->GetYAt(c->GetX(), bz, xyz[1]);
1125 par->GetZAt(c->GetX(), bz, xyz[2]);
1126 vecX[i] = c->GetX();
1132 if (fStreamLevel>5 && count1>10 && count2>10){
1134 // huge output - cluster residuals to be investigated
1136 TTreeSRedirector *cstream = GetDebugStreamer();
1137 AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
1138 AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
1141 Track->Draw("Cl[].fY-vtY.fElements:vtY.fElements-vtX.fElements*tan(pi/18.)>>his(100,-10,0)","Cl.fY!=0&&abs(Cl.fY-vtY.fElements)<1","prof");
1146 (*cstream)<<"Track"<<
1147 "run="<<fRun<< // run number
1148 "event="<<fEvent<< // event number
1149 "time="<<fTime<< // time stamp of event
1150 "trigger="<<fTrigger<< // trigger
1151 "triggerClass="<<&fTriggerClass<< // trigger
1152 "mag="<<fMagF<< // magnetic field
1174 void AliTPCcalibAlign::Process12(const Double_t *t1,
1176 TLinearFitter *fitter) const
1178 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1179 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
1180 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
1181 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1182 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1184 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
1185 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
1186 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
1190 const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1191 const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1194 Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1195 Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1196 Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1197 Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1202 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1203 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
1204 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1205 for (Int_t i=0; i<12;i++) p[i]=0.;
1210 p[0+1] = y1*dydx2; // a01
1211 p[0+2] = z1*dydx2; // a02
1212 p[9+0] = dydx2; // a03
1214 fitter->AddPoint(p,value,sy);
1216 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1217 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
1218 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1219 for (Int_t i=0; i<12;i++) p[i]=0.;
1224 p[0+1] = y1*dzdx2; // a01
1225 p[0+2] = z1*dzdx2; // a02
1226 p[9+0] = dzdx2; // a03
1228 fitter->AddPoint(p,value,sz);
1230 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1231 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1232 for (Int_t i=0; i<12;i++) p[i]=0.;
1234 p[3+1] = dydx1; // a11
1235 p[3+2] = dzdx1; // a12
1236 p[0+0] = -dydx2; // a00
1237 p[0+1] = -dydx1*dydx2; // a01
1238 p[0+2] = -dzdx1*dydx2; // a02
1240 fitter->AddPoint(p,value,sdydx);
1242 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1243 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1244 for (Int_t i=0; i<12;i++) p[i]=0.;
1246 p[6+1] = dydx1; // a21
1247 p[6+2] = dzdx1; // a22
1248 p[0+0] = -dzdx2; // a00
1249 p[0+1] = -dydx1*dzdx2; // a01
1250 p[0+2] = -dzdx1*dzdx2; // a02
1252 fitter->AddPoint(p,value,sdzdx);
1255 void AliTPCcalibAlign::Process9(const Double_t * const t1,
1256 const Double_t * const t2,
1257 TLinearFitter *fitter) const
1259 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1260 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
1261 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
1262 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1263 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1265 // a00 a01 a02 a03 1 p[0] p[1] p[6]
1266 // a10 a11 a12 a13 ==> p[2] 1 p[3] p[7]
1267 // a20 a21 a21 a23 p[4] p[5] 1 p[8]
1270 const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1271 const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1273 Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1274 Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1275 Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1276 Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1282 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1283 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
1284 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1285 for (Int_t i=0; i<12;i++) p[i]=0.;
1290 p[0] += y1*dydx2; // a01
1291 p[1] += z1*dydx2; // a02
1292 p[6] += dydx2; // a03
1293 value = y2-y1; //-a11
1294 fitter->AddPoint(p,value,sy);
1296 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1297 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
1298 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1299 for (Int_t i=0; i<12;i++) p[i]=0.;
1304 p[0] += y1*dzdx2; // a01
1305 p[1] += z1*dzdx2; // a02
1306 p[6] += dzdx2; // a03
1307 value = z2-z1; //-a22
1308 fitter->AddPoint(p,value,sz);
1310 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1311 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1312 for (Int_t i=0; i<12;i++) p[i]=0.;
1314 //p[] += dydx1; // a11
1315 p[3] += dzdx1; // a12
1316 //p[] += -dydx2; // a00
1317 p[0] += -dydx1*dydx2; // a01
1318 p[1] += -dzdx1*dydx2; // a02
1319 value = -dydx1+dydx2; // -a11 + a00
1320 fitter->AddPoint(p,value,sdydx);
1322 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1323 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1324 for (Int_t i=0; i<12;i++) p[i]=0.;
1326 p[5] += dydx1; // a21
1327 //p[] += dzdx1; // a22
1328 //p[] += -dzdx2; // a00
1329 p[0] += -dydx1*dzdx2; // a01
1330 p[1] += -dzdx1*dzdx2; // a02
1331 value = -dzdx1+dzdx2; // -a22 + a00
1332 fitter->AddPoint(p,value,sdzdx);
1335 void AliTPCcalibAlign::Process6(const Double_t *const t1,
1336 const Double_t *const t2,
1337 TLinearFitter *fitter) const
1339 // x2 = 1 *x1 +-a01*y1 + 0 +a03
1340 // y2 = a01*x1 + 1 *y1 + 0 +a13
1341 // z2 = a20*x1 + a21*y1 + 1 *z1 +a23
1342 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1343 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
1345 // a00 a01 a02 a03 1 -p[0] 0 p[3]
1346 // a10 a11 a12 a13 ==> p[0] 1 0 p[4]
1347 // a20 a21 a21 a23 p[1] p[2] 1 p[5]
1349 const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
1350 const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
1353 Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
1354 Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
1355 Double_t sz = TMath::Sqrt(t1[8]*t1[8]+t2[8]*t2[8]);
1356 Double_t sdzdx = TMath::Sqrt(t1[9]*t1[9]+t2[9]*t2[9]);
1361 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1362 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
1363 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
1364 for (Int_t i=0; i<12;i++) p[i]=0.;
1369 p[0] += -y1*dydx2; // a01
1370 //p[] += z1*dydx2; // a02
1371 p[3] += dydx2; // a03
1372 value = y2-y1; //-a11
1373 fitter->AddPoint(p,value,sy);
1375 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
1376 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
1377 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
1378 for (Int_t i=0; i<12;i++) p[i]=0.;
1383 p[0] += -y1*dzdx2; // a01
1384 //p[] += z1*dzdx2; // a02
1385 p[3] += dzdx2; // a03
1386 value = z2-z1; //-a22
1387 fitter->AddPoint(p,value,sz);
1389 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1390 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
1391 for (Int_t i=0; i<12;i++) p[i]=0.;
1393 //p[] += dydx1; // a11
1394 //p[] += dzdx1; // a12
1395 //p[] += -dydx2; // a00
1396 //p[0] += dydx1*dydx2; // a01 FIXME- 0912 MI
1397 //p[] += -dzdx1*dydx2; // a02
1398 value = -dydx1+dydx2; // -a11 + a00
1399 fitter->AddPoint(p,value,sdydx);
1401 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
1402 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
1403 for (Int_t i=0; i<12;i++) p[i]=0.;
1405 // p[2] += dydx1; // a21 FIXME- 0912 MI
1406 //p[] += dzdx1; // a22
1407 //p[] += -dzdx2; // a00
1408 //p[0] += dydx1*dzdx2; // a01 FIXME- 0912 MI
1409 //p[] += -dzdx1*dzdx2; // a02
1410 value = -dzdx1+dzdx2; // -a22 + a00
1411 fitter->AddPoint(p,value,sdzdx);
1417 void AliTPCcalibAlign::EvalFitters(Int_t minPoints) {
1421 // Perform the fitting using linear fitters
1424 TFile fff("alignDebug.root","recreate");
1425 for (Int_t s1=0;s1<72;++s1)
1426 for (Int_t s2=0;s2<72;++s2){
1427 if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1428 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1430 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1431 f->Write(Form("f12_%d_%d",s1,s2));
1433 f->Write(Form("f12_%d_%d",s1,s2));
1436 if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1437 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1439 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1441 f->Write(Form("f9_%d_%d",s1,s2));
1444 if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>minPoints) {
1445 // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
1447 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
1449 f->Write(Form("f6_%d_%d",s1,s2));
1454 for (Int_t s1=0;s1<72;++s1)
1455 for (Int_t s2=0;s2<72;++s2){
1456 if (GetTransformation12(s1,s2,mat)){
1457 fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2));
1459 if (GetTransformation9(s1,s2,mat)){
1460 fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2));
1462 if (GetTransformation6(s1,s2,mat)){
1463 fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2));
1466 //this->Write("align");
1470 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
1472 // get or make fitter - general linear transformation
1474 static Int_t counter12=0;
1475 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]");
1476 TLinearFitter * fitter = GetFitter12(s1,s2);
1477 if (fitter) return fitter;
1478 // 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]");
1479 fitter =new TLinearFitter(&f12,"");
1480 fitter->StoreData(kFALSE);
1481 fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
1483 // if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter12<<endl;
1487 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
1489 //get or make fitter - general linear transformation - no scaling
1491 static Int_t counter9=0;
1492 static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
1493 TLinearFitter * fitter = GetFitter9(s1,s2);
1494 if (fitter) return fitter;
1495 // fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
1496 fitter =new TLinearFitter(&f9,"");
1497 fitter->StoreData(kFALSE);
1498 fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
1500 // if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter9<<endl;
1504 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
1506 // get or make fitter - 6 paramater linear tranformation
1509 // - tilting x-z, y-z
1510 static Int_t counter6=0;
1511 static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1512 TLinearFitter * fitter = GetFitter6(s1,s2);
1513 if (fitter) return fitter;
1514 // fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
1515 fitter=new TLinearFitter(&f6,"");
1516 fitter->StoreData(kFALSE);
1517 fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
1519 // if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<" : "<<counter6<<endl;
1527 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
1529 // GetTransformation matrix - 12 paramaters - generael linear transformation
1531 if (!GetFitter12(s1,s2))
1535 GetFitter12(s1,s2)->GetParameters(p);
1537 a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
1538 a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
1539 a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
1540 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
1545 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
1547 // GetTransformation matrix - 9 paramaters - general linear transformation
1550 if (!GetFitter9(s1,s2))
1554 GetFitter9(s1,s2)->GetParameters(p);
1556 a[0][0]=1; a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
1557 a[1][0]=p[2]; a[1][1]=1; a[1][2]=p[3]; a[1][3]=p[7];
1558 a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1; a[2][3]=p[8];
1559 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
1564 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
1566 // GetTransformation matrix - 6 paramaters
1569 // 2 tilting x-z y-z
1570 if (!GetFitter6(s1,s2))
1574 GetFitter6(s1,s2)->GetParameters(p);
1576 a[0][0]=1; a[0][1]=-p[0];a[0][2]=0; a[0][3]=p[3];
1577 a[1][0]=p[0]; a[1][1]=1; a[1][2]=0; a[1][3]=p[4];
1578 a[2][0]=p[1]; a[2][1]=p[2]; a[2][2]=1; a[2][3]=p[5];
1579 a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
1584 void AliTPCcalibAlign::MakeResidualHistos(){
1586 // Make cluster residual histograms
1588 Double_t xminTrack[9], xmaxTrack[9];
1590 TString axisName[9],axisTitle[9];
1592 // 0 - delta of interest
1593 // 1 - global phi in sector number as float
1598 axisName[0]="delta"; axisTitle[0]="#Delta (cm)";
1599 if (TMath::Abs(AliTracker::GetBz())<0.01){
1600 binsTrack[0]=60; xminTrack[0]=-1.2; xmaxTrack[0]=1.2;
1602 binsTrack[0]=60; xminTrack[0]=-0.6; xmaxTrack[0]=0.6;
1605 axisName[1]="sector"; axisTitle[1]="Sector Number";
1606 binsTrack[1]=180; xminTrack[1]=0; xmaxTrack[1]=18;
1608 axisName[2]="R"; axisTitle[2]="r (cm)";
1609 binsTrack[2]=53; xminTrack[2]=85.; xmaxTrack[2]=245.;
1612 axisName[3]="kZ"; axisTitle[3]="dz/dx";
1613 binsTrack[3]=36; xminTrack[3]=-1.8; xmaxTrack[3]=1.8;
1615 fClusterDelta[0] = new THnF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1616 fClusterDelta[1] = new THnF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1620 for (Int_t ivar=0;ivar<2;ivar++){
1621 for (Int_t ivar2=0;ivar2<4;ivar2++){
1622 fClusterDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1623 fClusterDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
1630 void AliTPCcalibAlign::MakeResidualHistosTracklet(){
1632 // Make tracklet residual histograms
1634 Double_t xminTrack[9], xmaxTrack[9];
1636 TString axisName[9],axisTitle[9];
1638 // 0 - delta of interest
1639 // 1 - global phi in sector number as float
1647 axisName[0]="delta"; axisTitle[0]="#Delta (cm)";
1648 binsTrack[0]=60; xminTrack[0]=-0.5; xmaxTrack[0]=0.5;
1650 axisName[1]="phi"; axisTitle[1]="#phi";
1651 binsTrack[1]=90; xminTrack[1]=-TMath::Pi(); xmaxTrack[1]=TMath::Pi();
1653 axisName[2]="localX"; axisTitle[2]="x (cm)";
1654 binsTrack[2]=10; xminTrack[2]=120.; xmaxTrack[2]=200.;
1656 axisName[3]="kY"; axisTitle[3]="dy/dx";
1657 binsTrack[3]=10; xminTrack[3]=-0.5; xmaxTrack[3]=0.5;
1659 axisName[4]="kZ"; axisTitle[4]="dz/dx";
1660 binsTrack[4]=11; xminTrack[4]=-1.1; xmaxTrack[4]=1.1;
1662 axisName[5]="is1"; axisTitle[5]="is1";
1663 binsTrack[5]=72; xminTrack[5]=0; xmaxTrack[5]=72;
1665 axisName[6]="is0"; axisTitle[6]="is0";
1666 binsTrack[6]=72; xminTrack[6]=0; xmaxTrack[6]=72;
1668 axisName[7]="z"; axisTitle[7]="z(cm)";
1669 binsTrack[7]=12; xminTrack[7]=-240; xmaxTrack[7]=240;
1671 axisName[8]="IsPrimary"; axisTitle[8]="Is Primary";
1672 binsTrack[8]=2; xminTrack[8]=-0.1; xmaxTrack[8]=1.1;
1675 xminTrack[0]=-0.25; xmaxTrack[0]=0.25;
1676 fTrackletDelta[0] = new THnSparseF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
1677 xminTrack[0]=-0.5; xmaxTrack[0]=0.5;
1678 fTrackletDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
1679 xminTrack[0]=-0.005; xmaxTrack[0]=0.005;
1680 fTrackletDelta[2] = new THnSparseF("#Delta_{kY}","#Delta_{kY}", 9, binsTrack,xminTrack, xmaxTrack);
1681 xminTrack[0]=-0.008; xmaxTrack[0]=0.008;
1682 fTrackletDelta[3] = new THnSparseF("#Delta_{kZ}","#Delta_{kZ}", 9, binsTrack,xminTrack, xmaxTrack);
1686 for (Int_t ivar=0;ivar<4;ivar++){
1687 for (Int_t ivar2=0;ivar2<9;ivar2++){
1688 fTrackletDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1689 fTrackletDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
1697 void AliTPCcalibAlign::FillHisto(const Double_t *t1,
1699 Int_t s1,Int_t s2) {
1701 // Fill residual histograms
1707 Double_t dy = t2[1]-t1[1];
1708 Double_t dphi = t2[2]-t1[2];
1709 Double_t dz = t2[3]-t1[3];
1710 Double_t dtheta = t2[4]-t1[4];
1711 Double_t zmean = (t2[3]+t1[3])*0.5;
1713 GetHisto(kPhi,s1,s2,kTRUE)->Fill(dphi);
1714 GetHisto(kTheta,s1,s2,kTRUE)->Fill(dtheta);
1715 GetHisto(kY,s1,s2,kTRUE)->Fill(dy);
1716 GetHisto(kZ,s1,s2,kTRUE)->Fill(dz);
1718 GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(zmean,dphi);
1719 GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(zmean,dtheta);
1720 GetHisto(kYz,s1,s2,kTRUE)->Fill(zmean,dy);
1721 GetHisto(kZz,s1,s2,kTRUE)->Fill(zmean,dz);
1723 GetHisto(kYPhi,s1,s2,kTRUE)->Fill(t2[2],dy);
1724 GetHisto(kZTheta,s1,s2,kTRUE)->Fill(t2[4],dz);
1729 void AliTPCcalibAlign::FillHisto(AliExternalTrackParam *tp1,
1730 AliExternalTrackParam *tp2,
1731 Int_t s1,Int_t s2) {
1733 // Fill residual histograms
1735 if (s2<s1) return;//
1736 const Double_t kEpsilon=0.001;
1737 Double_t x[9]={0,0,0,0,0,0,0,0,0};
1738 AliExternalTrackParam p1(*tp1);
1739 AliExternalTrackParam p2(*tp2);
1741 // inner outer - match at the IROC-OROC boundary
1742 if (!p1.PropagateTo(fXIO, AliTrackerBase::GetBz())) return;
1744 if (!p2.Rotate(p1.GetAlpha())) return;
1745 if (!p2.PropagateTo(p1.GetX(),AliTrackerBase::GetBz())) return;
1746 if (TMath::Abs(p1.GetX()-p2.GetX())>kEpsilon) return;
1749 x[1]=TMath::ATan2(xyz[1],xyz[0]);
1751 x[3]=0.5*(p1.GetSnp()+p2.GetSnp()); // mean snp
1752 x[4]=0.5*(p1.GetTgl()+p2.GetTgl()); // mean tgl
1755 x[7]=0.5*(p1.GetZ()+p2.GetZ());
1757 Int_t isPrimary = (TMath::Abs(p1.GetTgl()-p1.GetZ()/p1.GetX())<0.1) ? 1:0;
1760 x[0]=p2.GetY()-p1.GetY();
1761 fTrackletDelta[0]->Fill(x);
1762 x[0]=p2.GetZ()-p1.GetZ();
1763 fTrackletDelta[1]->Fill(x);
1764 x[0]=p2.GetSnp()-p1.GetSnp();
1765 fTrackletDelta[2]->Fill(x);
1766 x[0]=p2.GetTgl()-p1.GetTgl();
1767 fTrackletDelta[3]->Fill(x);
1768 TTreeSRedirector *cstream = GetDebugStreamer();
1770 (*cstream)<<"trackletMatch"<<
1771 "tp1.="<<tp1<< // input tracklet
1773 "p1.="<<&p1<< // tracklet in the ref frame
1784 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
1787 // return specified residual histogram - it is only QA
1788 // if force specified the histogram and given histogram is not existing
1789 // new histogram is created
1791 if (GetIndex(s1,s2)>=72*72) return 0;
1792 TObjArray *histoArray=0;
1795 histoArray = &fDyHistArray; break;
1797 histoArray = &fDzHistArray; break;
1799 histoArray = &fDphiHistArray; break;
1801 histoArray = &fDthetaHistArray; break;
1803 histoArray = &fDyPhiHistArray; break;
1805 histoArray = &fDzThetaHistArray; break;
1807 histoArray = &fDyZHistArray; break;
1809 histoArray = &fDzZHistArray; break;
1811 histoArray = &fDphiZHistArray; break;
1813 histoArray = &fDthetaZHistArray; break;
1815 TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
1816 if (histo) return histo;
1817 if (force==kFALSE) return 0;
1823 name<<"hist_y_"<<s1<<"_"<<s2;
1824 title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1825 histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.5,0.5); // +/- 5 mm
1828 name<<"hist_z_"<<s1<<"_"<<s2;
1829 title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1830 histo = new TH1D(name.str().c_str(),title.str().c_str(),100,-0.3,0.3); // +/- 3 mm
1833 name<<"hist_phi_"<<s1<<"_"<<s2;
1834 title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1835 histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.01,0.01); // +/- 10 mrad
1838 name<<"hist_theta_"<<s1<<"_"<<s2;
1839 title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1840 histo =new TH1D(name.str().c_str(),title.str().c_str(),100,-0.01,0.01); // +/- 10 mrad
1845 name<<"hist_yphi_"<<s1<<"_"<<s2;
1846 title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1847 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,100,-0.5,0.5); // +/- 5 mm
1850 name<<"hist_ztheta_"<<s1<<"_"<<s2;
1851 title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1852 histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,100,-0.3,0.3); // +/- 3 mm
1858 name<<"hist_yz_"<<s1<<"_"<<s2;
1859 title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1860 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.5,0.5); // +/- 5 mm
1863 name<<"hist_zz_"<<s1<<"_"<<s2;
1864 title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1865 histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.3,0.3); // +/- 3 mm
1868 name<<"hist_phiz_"<<s1<<"_"<<s2;
1869 title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1870 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.01,0.01); // +/- 10 mrad
1873 name<<"hist_thetaz_"<<s1<<"_"<<s2;
1874 title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1875 histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,100,-0.01,0.01); // +/- 10 mrad
1880 histo->SetDirectory(0);
1881 histoArray->AddAt(histo,GetIndex(s1,s2));
1885 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec,
1886 Int_t i0, Int_t i1, FitType type)
1892 //TObjArray *fitArray=0;
1893 Double_t xsec[1000];
1894 Double_t ysec[1000];
1896 for (Int_t isec = sec0; isec<=sec1; isec++){
1897 Int_t isec2 = (isec+dsec)%72;
1900 GetTransformation6(isec,isec2,mat);break;
1902 GetTransformation9(isec,isec2,mat);break;
1904 GetTransformation12(isec,isec2,mat);break;
1907 ysec[npoints]=mat(i0,i1);
1910 TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1912 snprintf(name,100,"Mat[%d,%d] Type=%d",i0,i1,type);
1917 void AliTPCcalibAlign::MakeTree(const char *fname, Int_t minPoints){
1919 // make tree with alignment cosntant -
1920 // For QA visualization
1923 TFile fcalib("CalibObjects.root");
1924 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
1925 AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
1926 align->EvalFitters();
1927 align->MakeTree("alignTree.root");
1928 TFile falignTree("alignTree.root");
1929 TTree * treeAlign = (TTree*)falignTree.Get("Align");
1931 TTreeSRedirector cstream(fname);
1932 for (Int_t s1=0;s1<72;++s1)
1933 for (Int_t s2=0;s2<72;++s2){
1938 TVectorD param6Diff; // align parameters diff
1939 TVectorD param6s1(6); // align parameters sector1
1940 TVectorD param6s2(6); // align parameters sector2
1945 TMatrixD * kpar = fSectorParamA;
1946 TMatrixD * kcov = fSectorCovarA;
1948 kpar = fSectorParamC;
1949 kcov = fSectorCovarC;
1951 for (Int_t ipar=0;ipar<6;ipar++){
1952 Int_t isec1 = s1%18;
1953 Int_t isec2 = s2%18;
1954 if (s1>35) isec1+=18;
1955 if (s2>35) isec2+=18;
1956 param6s1(ipar)=(*kpar)(6*isec1+ipar,0);
1957 param6s2(ipar)=(*kpar)(6*isec2+ipar,0);
1961 Double_t dy=0, dz=0, dphi=0,dtheta=0;
1962 Double_t sy=0, sz=0, sphi=0,stheta=0;
1963 Double_t ny=0, nz=0, nphi=0,ntheta=0;
1964 Double_t chi2v12=0, chi2v9=0, chi2v6=0;
1966 // TLinearFitter * fitter = 0;
1967 if (fPoints[GetIndex(s1,s2)]>minPoints){
1971 // fitter = GetFitter12(s1,s2);
1972 // npoints = fitter->GetNpoints();
1973 // chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1976 // fitter = GetFitter9(s1,s2);
1977 // npoints = fitter->GetNpoints();
1978 // chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1980 // fitter = GetFitter6(s1,s2);
1981 // npoints = fitter->GetNpoints();
1982 // chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1983 // fitter->GetParameters(param6Diff);
1985 // GetTransformation6(s1,s2,m6);
1986 // GetTransformation9(s1,s2,m9);
1987 // GetTransformation12(s1,s2,m12);
1989 // fitter = GetFitter6(s1,s2);
1990 // //fitter->FixParameter(3,0);
1991 // //fitter->Eval();
1992 // GetTransformation6(s1,s2,m6FX);
1995 his = GetHisto(kY,s1,s2);
1996 if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1997 his = GetHisto(kZ,s1,s2);
1998 if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1999 his = GetHisto(kPhi,s1,s2);
2000 if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
2001 his = GetHisto(kTheta,s1,s2);
2002 if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
2006 AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
2007 if (!magF) AliError("Magneticd field - not initialized");
2008 Double_t bz = magF->SolenoidField()/10.; //field in T
2011 "run="<<fRun<< // run
2013 "s1="<<s1<< // reference sector
2014 "s2="<<s2<< // sector to align
2015 "m6FX.="<<&m6FX<< // tranformation matrix
2016 "m6.="<<&m6<< // tranformation matrix
2019 "chi2v12="<<chi2v12<<
2023 "p6.="<<¶m6Diff<<
2024 "p6s1.="<<¶m6s1<<
2025 "p6s2.="<<¶m6s2<<
2026 // histograms mean RMS and entries
2045 //_____________________________________________________________________
2046 Long64_t AliTPCcalibAlign::Merge(TCollection* const list) {
2050 if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
2053 if (list->IsEmpty())
2056 TIterator* iter = list->MakeIterator();
2061 TString str1(GetName());
2062 while((obj = iter->Next()) != 0)
2064 AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
2065 if (entry == 0) continue;
2066 if (str1.CompareTo(entry->GetName())!=0) continue;
2074 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
2076 // Add entry - used for merging of compoents
2078 for (Int_t i=0; i<72;i++){
2079 for (Int_t j=0; j<72;j++){
2080 if (align->fPoints[GetIndex(i,j)]<1) continue;
2081 fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
2085 for (Int_t itype=0; itype<10; itype++){
2086 TH1 * his0=0, *his1=0;
2087 his0 = GetHisto((HistoType)itype,i,j);
2088 his1 = align->GetHisto((HistoType)itype,i,j);
2090 if (his0) his0->Add(his1);
2092 his0 = GetHisto((HistoType)itype,i,j,kTRUE);
2099 TLinearFitter *f0=0;
2100 TLinearFitter *f1=0;
2101 for (Int_t i=0; i<72;i++){
2102 for (Int_t j=0; j<72;j++){
2103 if (align->fPoints[GetIndex(i,j)]<1) continue;
2107 f0 = GetFitter12(i,j);
2108 f1 = align->GetFitter12(i,j);
2110 if (f0) f0->Add(f1);
2112 fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
2117 f0 = GetFitter9(i,j);
2118 f1 = align->GetFitter9(i,j);
2120 if (f0) f0->Add(f1);
2122 fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
2125 f0 = GetFitter6(i,j);
2126 f1 = align->GetFitter6(i,j);
2128 if (f0) f0->Add(f1);
2130 fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
2136 // Add Kalman filter
2138 for (Int_t i=0;i<36;i++){
2139 TMatrixD *par0 = (TMatrixD*)fArraySectorIntParam.At(i);
2142 par0 = (TMatrixD*)fArraySectorIntParam.At(i);
2144 TMatrixD *par1 = (TMatrixD*)align->fArraySectorIntParam.At(i);
2145 if (!par1) continue;
2147 TMatrixD *cov0 = (TMatrixD*)fArraySectorIntCovar.At(i);
2148 TMatrixD *cov1 = (TMatrixD*)align->fArraySectorIntCovar.At(i);
2149 UpdateSectorKalman(*par0,*cov0,*par1,*cov1);
2151 if (!fSectorParamA){
2154 if (align->fSectorParamA){
2155 UpdateKalman(*fSectorParamA,*fSectorCovarA,*align->fSectorParamA,*align->fSectorCovarA);
2156 UpdateKalman(*fSectorParamC,*fSectorCovarC,*align->fSectorParamC,*align->fSectorCovarC);
2158 if (!fClusterDelta[0]) MakeResidualHistos();
2160 for (Int_t i=0; i<2; i++){
2161 if (align->fClusterDelta[i]){
2162 fClusterDelta[i]->Add(align->fClusterDelta[i]);
2167 for (Int_t i=0; i<4; i++){
2168 if (!fTrackletDelta[i] && align->fTrackletDelta[i]) {
2169 fTrackletDelta[i]= (THnSparse*)(align->fTrackletDelta[i]->Clone());
2172 if (align->fTrackletDelta[i]) {
2173 if (fTrackletDelta[i]->GetEntries()<fgkMergeEntriesCut){
2174 fTrackletDelta[i]->Add(align->fTrackletDelta[i]);
2181 Double_t AliTPCcalibAlign::Correct(Int_t type, Int_t value, Int_t s1, Int_t s2, Double_t x1, Double_t y1, Double_t z1, Double_t dydx1,Double_t dzdx1){
2183 // GetTransformed value
2186 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
2187 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
2188 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
2189 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
2190 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
2193 const TMatrixD * mat = GetTransformation(s1,s2,type);
2195 if (value==0) return x1;
2196 if (value==1) return y1;
2197 if (value==2) return z1;
2198 if (value==3) return dydx1;
2199 if (value==4) return dzdx1;
2201 if (value==5) return dydx1;
2202 if (value==6) return dzdx1;
2208 valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
2212 valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
2215 valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
2218 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
2219 valT = (*mat)(1,0) +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
2220 valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
2224 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
2225 valT = (*mat)(2,0) +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
2226 valT/= ((*mat)(0,0) +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
2230 // onlys shift in angle
2231 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
2232 valT = (*mat)(1,0) +(*mat)(1,1)*dydx1;
2236 // only shift in angle
2237 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
2238 valT = (*mat)(2,0) +(*mat)(2,1)*dydx1;
2245 void AliTPCcalibAlign::Constrain1Pt(AliExternalTrackParam &track1, const AliExternalTrackParam &track2, Bool_t noField){
2247 // Update track parameters t1
2249 TMatrixD vecXk(5,1); // X vector
2250 TMatrixD covXk(5,5); // X covariance
2251 TMatrixD matHk(1,5); // vector to mesurement
2252 TMatrixD measR(1,1); // measurement error
2253 //TMatrixD matQk(5,5); // prediction noise vector
2254 TMatrixD vecZk(1,1); // measurement
2256 TMatrixD vecYk(1,1); // Innovation or measurement residual
2257 TMatrixD matHkT(5,1);
2258 TMatrixD matSk(1,1); // Innovation (or residual) covariance
2259 TMatrixD matKk(5,1); // Optimal Kalman gain
2260 TMatrixD mat1(5,5); // update covariance matrix
2261 TMatrixD covXk2(5,5); //
2262 TMatrixD covOut(5,5);
2264 Double_t *param1=(Double_t*) track1.GetParameter();
2265 Double_t *covar1=(Double_t*) track1.GetCovariance();
2268 // copy data to the matrix
2269 for (Int_t ipar=0; ipar<5; ipar++){
2270 vecXk(ipar,0) = param1[ipar];
2271 for (Int_t jpar=0; jpar<5; jpar++){
2272 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
2278 vecZk(0,0) = track2.GetParameter()[4]; // 1/pt measurement from track 2
2279 measR(0,0) = track2.GetCovariance()[14]; // 1/pt measurement error
2281 measR(0,0)*=0.000000001;
2285 matHk(0,0)=0; matHk(0,1)= 0; matHk(0,2)= 0;
2286 matHk(0,3)= 0; matHk(0,4)= 1; // vector to measurement
2290 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
2291 matHkT=matHk.T(); matHk.T();
2292 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
2294 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
2295 vecXk += matKk*vecYk; // updated vector
2296 mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
2297 covXk2 = (mat1-(matKk*matHk));
2298 covOut = covXk2*covXk;
2302 // copy from matrix to parameters
2311 for (Int_t ipar=0; ipar<5; ipar++){
2312 param1[ipar]= vecXk(ipar,0) ;
2313 for (Int_t jpar=0; jpar<5; jpar++){
2314 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
2320 void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t niter){
2322 // Global Align -combine the partial alignment of pair of sectors
2323 // minPoints - minimal number of points - don't use sector alignment wit smaller number
2324 // sysError - error added to the alignemnt error
2326 AliTPCcalibAlign * align = this;
2327 TMatrixD * arrayAlign[72];
2328 TMatrixD * arrayAlignDiff[72];
2330 for (Int_t i=0;i<72; i++) {
2331 TMatrixD * mat = new TMatrixD(4,4);
2334 arrayAlignDiff[i]=(TMatrixD*)(mat->Clone());
2337 TTreeSRedirector *cstream = new TTreeSRedirector("galign6.root");
2338 for (Int_t iter=0; iter<niter;iter++){
2339 printf("Iter=\t%d\n",iter);
2340 for (Int_t is0=0;is0<72; is0++) {
2342 //TMatrixD *mati0 = arrayAlign[is0];
2343 TMatrixD matDiff(4,4);
2345 for (Int_t is1=0;is1<72; is1++) {
2346 Bool_t invers=kFALSE;
2350 const TMatrixD *mat = align->GetTransformation(is0,is1,0);
2352 npoints = align->GetFitter6(is0,is1)->GetNpoints();
2353 if (npoints>minPoints){
2354 align->GetFitter6(is0,is1)->GetCovarianceMatrix(covar);
2355 align->GetFitter6(is0,is1)->GetErrors(errors);
2360 mat = align->GetTransformation(is1,is0,0);
2362 npoints = align->GetFitter6(is1,is0)->GetNpoints();
2363 if (npoints>minPoints){
2364 align->GetFitter6(is1,is0)->GetCovarianceMatrix(covar);
2365 align->GetFitter6(is1,is0)->GetErrors(errors);
2370 if (npoints<minPoints) continue;
2373 if (is1/36>is0/36) weight*=2./3.; //IROC-OROC
2374 if (is1/36<is0/36) weight*=1./3.; //OROC-IROC
2375 if (is1/36==is0/36) weight*=1/3.; //OROC-OROC
2376 if (is1%36!=is0%36) weight*=1/2.; //Not up-down
2377 weight/=(errors[4]*errors[4]+sysError*sysError); // wieghting with error in Y
2380 TMatrixD matT = *mat;
2381 if (invers) matT.Invert();
2382 TMatrixD diffMat= (*(arrayAlign[is1]))*matT;
2383 diffMat-=(*arrayAlign[is0]);
2384 matDiff+=weight*diffMat;
2387 (*cstream)<<"LAlign"<<
2391 "npoints="<<npoints<<
2392 "m60.="<<arrayAlign[is0]<<
2393 "m61.="<<arrayAlign[is1]<<
2395 "diff.="<<&diffMat<<
2406 (*arrayAlignDiff[is0]) = matDiff;
2409 for (Int_t is0=0;is0<72; is0++) {
2410 if (is0<36) (*arrayAlign[is0]) += 0.4*(*arrayAlignDiff[is0]);
2411 if (is0>=36) (*arrayAlign[is0]) += 0.2*(*arrayAlignDiff[is0]);
2413 (*cstream)<<"GAlign"<<
2416 "m6.="<<arrayAlign[is0]<<
2422 for (Int_t isec=0;isec<72;isec++){
2423 fCombinedMatrixArray6.AddAt(arrayAlign[isec],isec);
2424 delete arrayAlignDiff[isec];
2429 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){
2431 // Refit tracklet linearly using clusters at given sector isec
2432 // Clusters are rotated to the reference frame of sector refSector
2434 // fit parameters and errors retruning in the fitParam
2436 // seed - acces to the original clusters
2437 // isec - sector to be refited
2449 // ref sector is the sector defining ref frame - rotation
2450 // return value - number of used clusters
2452 const Int_t kMinClusterF=15;
2453 const Int_t kdrow1 =10; // rows to skip at the end
2454 const Int_t kdrow0 =3; // rows to skip at beginning
2455 const Float_t kedgeyIn=2.5;
2456 const Float_t kedgeyOut=4.0;
2457 const Float_t kMaxDist=5; // max distance -in sigma
2458 const Float_t kMaxCorrY=0.05; // max correction
2460 Double_t dalpha = 0;
2461 if ((refSector%18)!=(isec%18)){
2462 dalpha = -((refSector%18)-(isec%18))*TMath::TwoPi()/18.;
2464 Double_t ca = TMath::Cos(dalpha);
2465 Double_t sa = TMath::Sin(dalpha);
2468 AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
2470 // full track fit parameters
2472 static TLinearFitter fyf(2,"pol1"); // change to static - suggestion of calgrind - 30 % of time
2473 static TLinearFitter fzf(2,"pol1"); // relative to time of given class
2474 TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
2475 TMatrixD covY(4,4),covZ(4,4);
2476 Double_t chi2FacY =1;
2477 Double_t chi2FacZ =1;
2482 Float_t erry=0.1; // initial cluster error estimate
2483 Float_t errz=0.1; // initial cluster error estimate
2484 for (Int_t iter=0; iter<2; iter++){
2487 for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
2488 AliTPCclusterMI *c=track->GetClusterPointer(irow);
2491 if (c->GetDetector()%36!=(isec%36)) continue;
2492 if (!both && c->GetDetector()!=isec) continue;
2494 if (c->GetRow()<kdrow0) continue;
2495 //cluster position in reference frame
2496 Double_t lxR = ca*c->GetX()-sa*c->GetY();
2497 Double_t lyR = +sa*c->GetX()+ca*c->GetY();
2498 Double_t lzR = c->GetZ();
2500 Double_t dx = lxR -xRef; // distance to reference X
2501 Double_t x[2]={dx, dx*dx};
2503 Double_t yfitR = pyf[0]+pyf[1]*dx; // fit value Y in ref frame
2504 Double_t zfitR = pzf[0]+pzf[1]*dx; // fit value Z in ref frame
2506 Double_t yfit = -sa*lxR + ca*yfitR; // fit value Y in local frame
2508 if (iter==0 &&c->GetType()<0) continue;
2510 if (TMath::Abs(lyR-yfitR)>kMaxDist*erry) continue;
2511 if (TMath::Abs(lzR-zfitR)>kMaxDist*errz) continue;
2512 Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2513 if (isec<36 && dedge<kedgeyIn) continue;
2514 if (isec>35 && dedge<kedgeyOut) continue;
2516 corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
2517 c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2519 corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
2520 c->GetY(),c->GetY(), c->GetZ(), pyf[1], c->GetMax(),2.5);
2521 if (TMath::Abs((corrtrY+corrclY)*0.5)>kMaxCorrY) continue;
2522 if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2524 fyf.AddPoint(x,lyR,erry);
2525 fzf.AddPoint(x,lzR,errz);
2527 nf = fyf.GetNpoints();
2528 if (nf<kMinClusterF) return 0; // not enough points - skip
2530 fyf.GetParameters(pyf);
2531 fyf.GetErrors(peyf);
2533 fzf.GetParameters(pzf);
2534 fzf.GetErrors(pezf);
2535 chi2FacY = TMath::Sqrt(fyf.GetChisquare()/(fyf.GetNpoints()-2.));
2536 chi2FacZ = TMath::Sqrt(fzf.GetChisquare()/(fzf.GetNpoints()-2.));
2543 fyf.GetCovarianceMatrix(covY);
2544 fzf.GetCovarianceMatrix(covZ);
2545 for (Int_t i0=0;i0<2;i0++)
2546 for (Int_t i1=0;i1<2;i1++){
2547 covY(i0,i1)*=chi2FacY*chi2FacY;
2548 covZ(i0,i1)*=chi2FacZ*chi2FacZ;
2553 fitParam[1] = pyf[0];
2554 fitParam[2] = pyf[1];
2555 fitParam[3] = pzf[0];
2556 fitParam[4] = pzf[1];
2559 fitParam[6] = peyf[0];
2560 fitParam[7] = peyf[1];
2561 fitParam[8] = pezf[0];
2562 fitParam[9] = pezf[1];
2565 tparam(0,0) = pyf[0];
2566 tparam(1,0) = pyf[1];
2567 tparam(2,0) = pzf[0];
2568 tparam(3,0) = pzf[1];
2570 tcovar(0,0) = covY(0,0);
2571 tcovar(1,1) = covY(1,1);
2572 tcovar(1,0) = covY(1,0);
2573 tcovar(0,1) = covY(0,1);
2574 tcovar(2,2) = covZ(0,0);
2575 tcovar(3,3) = covZ(1,1);
2576 tcovar(3,2) = covZ(1,0);
2577 tcovar(2,3) = covZ(0,1);
2581 void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){
2583 // Update the cluster residula histograms for setup with field
2584 // Kalman track fitting is used
2585 // Only high momenta primary tracks used
2587 // 1. Apply selection
2588 // 2. Refit the track - in-out
2589 // 3. Refit the track - out-in
2590 // 4. Combine In and Out track - - fil cluster residuals
2592 if (!fCurrentFriendTrack) return;
2593 if (!fCurrentFriendTrack->GetTPCOut()) return;
2594 const Double_t kPtCut=1.0; // pt
2595 const Double_t kSnpCut=0.2; // snp cut
2596 const Double_t kNclCut=120; //
2597 const Double_t kVertexCut=1;
2598 const Double_t kMaxDist=0.5; // max distance between tracks and cluster
2599 const Double_t kEdgeCut = 2.5;
2600 const Double_t kDelta2=0.2*0.2; // initial increase in covar matrix
2601 const Double_t kSigma=0.3; // error increase towards edges of TPC
2602 const Double_t kSkipBoundary=7.5; // skip track updates in the boundary IFC,OFC, IO
2604 if (!fCurrentTrack) return;
2605 if (!fCurrentFriendTrack) return;
2606 Float_t vertexXY=0,vertexZ=0;
2607 fCurrentTrack->GetImpactParameters(vertexXY,vertexZ);
2608 if (TMath::Abs(vertexXY)>kVertexCut) return;
2609 if (TMath::Abs(vertexZ)>kVertexCut) return;
2610 if (TMath::Abs(seed->Pt())<kPtCut) return;
2611 if (seed->GetNumberOfClusters()<kNclCut) return;
2612 if (TMath::Abs(seed->GetSnp())>kSnpCut) return;
2613 if (!fClusterDelta[0]) MakeResidualHistos();
2615 AliExternalTrackParam fitIn[160];
2616 AliExternalTrackParam fitOut[160];
2617 AliTPCROC * roc = AliTPCROC::Instance();
2618 Double_t xmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
2619 Double_t xDiff = ( -roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
2620 Double_t xIFC = ( roc->GetPadRowRadii(0,0));
2621 Double_t xOFC = ( roc->GetPadRowRadii(36,roc->GetNRows(36)-1));
2626 AliExternalTrackParam trackIn = *(fCurrentTrack->GetInnerParam());
2627 AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut());
2628 trackIn.ResetCovariance(10);
2629 trackOut.ResetCovariance(10);
2630 Double_t *covarIn = (Double_t*)trackIn.GetCovariance();
2631 Double_t *covarOut = (Double_t*)trackOut.GetCovariance();
2632 covarIn[0]+=kDelta2; covarIn[2]+=kDelta2;
2633 covarIn[5]+=kDelta2/(100.*100.); covarIn[9]=kDelta2/(100.*100.);
2634 covarIn[14]+=kDelta2/(5.*5.);
2635 covarOut[0]+=kDelta2; covarOut[2]+=kDelta2;
2636 covarOut[5]+=kDelta2/(100.*100.); covarOut[9]=kDelta2/(100.*100.);
2637 covarOut[14]+=kDelta2/(5.*5.);
2639 static Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2642 for (Int_t irow=0; irow<160; irow++){
2643 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2645 if (cl->GetX()<80) continue;
2646 if (detector<0) detector=cl->GetDetector()%36;
2647 if (detector!=cl->GetDetector()%36) return; // cluster from different sectors
2651 if (ncl<kNclCut) return;
2652 Int_t nclIn=0,nclOut=0;
2655 // Refit out - store residual maps
2657 for (Int_t irow=0; irow<160; irow++){
2658 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2660 if (cl->GetX()<80) continue;
2661 if (detector<0) detector=cl->GetDetector()%36;
2662 Int_t sector = cl->GetDetector();
2663 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
2664 if (cl->GetDetector()%36!=detector) continue;
2665 if (TMath::Abs(dalpha)>0.01){
2666 if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
2668 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
2669 Double_t cov[3]={0.1,0.,0.1};
2670 Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackOut.GetY());
2671 Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
2674 cov[0]+=kSigma*dmiddle; // bigger error at boundary
2675 cov[0]+=kSigma*dmiddle; // bigger error at boundary
2676 cov[2]+=kSigma*dmiddle; // bigger error at boundary
2677 cov[2]+=kSigma*dmiddle; // bigger error at boundary
2678 cov[0]+=kSigma/dedge; // bigger error close to the boundary
2679 cov[2]+=kSigma/dedge; // bigger error close to the boundary
2682 if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
2683 if (TMath::Abs(dedge)<kEdgeCut) continue;
2685 Bool_t doUpdate=kTRUE;
2686 if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
2687 if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
2688 if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
2690 if (TMath::Abs(cl->GetY()-trackOut.GetY())<kMaxDist){
2692 if (doUpdate) trackOut.Update(&r[1],cov);
2694 fitOut[irow]=trackOut;
2698 // Refit In - store residual maps
2700 for (Int_t irow=159; irow>=0; irow--){
2701 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2703 if (cl->GetX()<80) continue;
2704 if (detector<0) detector=cl->GetDetector()%36;
2705 Int_t sector = cl->GetDetector();
2706 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
2707 if (cl->GetDetector()%36!=detector) continue;
2708 if (TMath::Abs(dalpha)>0.01){
2709 if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
2711 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
2712 Double_t cov[3]={0.1,0.,0.1};
2713 Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackIn.GetY());
2714 Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
2717 cov[0]+=kSigma*dmiddle; // bigger error at boundary
2718 cov[0]+=kSigma*dmiddle; // bigger error at boundary
2719 cov[2]+=kSigma*dmiddle; // bigger error at boundary
2720 cov[2]+=kSigma*dmiddle; // bigger error at boundary
2721 cov[0]+=kSigma/dedge; // bigger error close to the boundary
2722 cov[2]+=kSigma/dedge; // bigger error close to the boundary
2725 if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
2726 if (TMath::Abs(dedge)<kEdgeCut) continue;
2727 Bool_t doUpdate=kTRUE;
2728 if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
2729 if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
2730 if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
2731 if (TMath::Abs(cl->GetY()-trackIn.GetY())<kMaxDist){
2733 if (doUpdate) trackIn.Update(&r[1],cov);
2735 fitIn[irow]=trackIn;
2739 for (Int_t irow=159; irow>=0; irow--){
2741 // Update kalman - +- direction
2742 // Store cluster residuals
2743 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
2745 if (cl->GetX()<80) continue;
2746 if (detector<0) detector=cl->GetDetector()%36;
2747 if (cl->GetDetector()%36!=detector) continue;
2748 if (fitIn[irow].GetX()<80) continue;
2749 if (fitOut[irow].GetX()<80) continue;
2750 AliExternalTrackParam trackSmooth = fitIn[irow];
2751 AliTrackerBase::UpdateTrack(trackSmooth, fitOut[irow]);
2753 Double_t resVector[5];
2754 trackSmooth.GetXYZ(xyz);
2755 resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
2756 if (resVector[1]<0) resVector[1]+=18;
2757 resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
2758 resVector[3]= cl->GetZ()/resVector[2];
2760 resVector[0]= cl->GetY()-trackSmooth.GetY();
2761 fClusterDelta[0]->Fill(resVector);
2762 resVector[0]= cl->GetZ()-trackSmooth.GetZ();
2763 fClusterDelta[1]->Fill(resVector);
2769 void AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
2771 // Update Kalman filter of Alignment - only setup without filed
2772 // IROC - OROC quadrants
2774 if (TMath::Abs(AliTracker::GetBz())>0.5) return;
2775 if (!fClusterDelta[0]) MakeResidualHistos();
2776 // const Int_t kMinClusterF=40;
2777 const Int_t kMinClusterFit=10;
2778 const Int_t kMinClusterQ=10;
2780 const Int_t kdrow1Fit =5; // rows to skip from fit at the end
2781 const Int_t kdrow0Fit =10; // rows to skip from fit at beginning
2782 const Float_t kedgey=3.0;
2783 const Float_t kMaxDist=1;
2784 const Float_t kMaxCorrY=0.05;
2785 const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF
2786 isec = isec%36; // use the hardware numbering
2789 AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
2791 // full track fit parameters
2793 static TLinearFitter fyf(2,"pol1"); // make it static - too much time for comiling of formula
2794 static TLinearFitter fzf(2,"pol1"); // calgrind recomendation
2795 TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
2796 TVectorD pyfc(2),pzfc(2);
2799 // make full fit as reference
2801 for (Int_t iter=0; iter<2; iter++){
2804 for (Int_t irow=kdrow0Fit;irow<159-kdrow1Fit;irow++) {
2805 AliTPCclusterMI *c=track->GetClusterPointer(irow);
2807 if ((c->GetDetector()%36)!=isec) continue;
2808 if (c->GetRow()<kdrow0Fit) continue;
2809 Double_t dx = c->GetX()-fXmiddle;
2810 Double_t x[2]={dx, dx*dx};
2811 if (iter==0 &&c->GetType()<0) continue;
2813 Double_t yfit = pyf[0]+pyf[1]*dx;
2814 Double_t zfit = pzf[0]+pzf[1]*dx;
2815 Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2816 if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
2817 if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
2818 if (dedge<kedgey) continue;
2820 corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
2821 c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2822 if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2824 if (TMath::Abs(x[0])<10){
2825 fyf.AddPoint(x,c->GetY(),0.1); //use only middle rows+-10cm
2826 fzf.AddPoint(x,c->GetZ(),0.1);
2829 nf = fyf.GetNpoints();
2830 if (fyf.GetNpoints()<kMinClusterFit) return; // not enough points - skip
2831 if (fzf.GetNpoints()<kMinClusterFit) return; // not enough points - skip
2833 fyf.GetParameters(pyf);
2834 fyf.GetErrors(peyf);
2836 fzf.GetParameters(pzf);
2837 fzf.GetErrors(pezf);
2842 TVectorD vecX(160); // x vector
2843 TVectorD vecY(160); // residuals vector
2844 TVectorD vecZ(160); // residuals vector
2845 TVectorD vPosG(3); //vertex position
2846 TVectorD vPosL(3); // vertex position in the TPC local system
2847 TVectorF vImpact(2); //track impact parameter
2848 // Double_t tofSignal= fCurrentTrack->GetTOFsignal(); // tof signal
2849 TVectorF tpcPosG(3); // global position of track at the middle of fXmiddle
2850 Double_t lphi = TMath::ATan2(pyf[0],fXmiddle); // expected local phi angle - if vertex at 0
2851 Double_t gphi = 2.*TMath::Pi()*(isec%18+0.5)/18.+lphi; // expected global phi if vertex at 0
2852 Double_t ky = pyf[0]/fXmiddle;
2853 Double_t kyE =0, kzE=0; // ky and kz expected
2854 Double_t alpha =2.*TMath::Pi()*(isec%18+0.5)/18.;
2855 Double_t scos=TMath::Cos(alpha);
2856 Double_t ssin=TMath::Sin(alpha);
2857 const AliESDVertex* vertex = fCurrentEvent->GetPrimaryVertexTracks();
2858 vertex->GetXYZ(vPosG.GetMatrixArray());
2859 fCurrentTrack->GetImpactParameters(vImpact[0],vImpact[1]); // track impact parameters
2861 tpcPosG[0]= scos*fXmiddle-ssin*pyf[0];
2862 tpcPosG[1]=+ssin*fXmiddle+scos*pyf[0];
2864 vPosL[0]= scos*vPosG[0]+ssin*vPosG[1];
2865 vPosL[1]=-ssin*vPosG[0]+scos*vPosG[1];
2867 kyE = (pyf[0]-vPosL[1])/(fXmiddle-vPosL[0]);
2868 kzE = (pzf[0]-vPosL[2])/(fXmiddle-vPosL[0]);
2870 // get constrained parameters
2872 Double_t xvertex=vPosL[0]-fXmiddle;
2873 fyf.AddPoint(&xvertex,vPosL[1], 0.00001);
2874 fzf.AddPoint(&xvertex,vPosL[2], 2.);
2876 fyf.GetParameters(pyfc);
2878 fzf.GetParameters(pzfc);
2881 // Make Fitters and params for 5 fitters
2882 // 1-4 OROC quadrants
2885 static TLinearFitter *fittersY[5]={0,0,0,0,0}; // calgrind recomendation - fater to clear points
2886 static TLinearFitter *fittersZ[5]={0,0,0,0,0}; // than create the fitter
2887 if (fittersY[0]==0){
2888 for (Int_t i=0;i<5;i++) {
2889 fittersY[i] = new TLinearFitter(2,"pol1");
2890 fittersZ[i] = new TLinearFitter(2,"pol1");
2895 TVectorD paramsY[5];
2896 TVectorD errorsY[5];
2899 TVectorD paramsZ[5];
2900 TVectorD errorsZ[5];
2903 for (Int_t i=0;i<5;i++) {
2905 paramsY[i].ResizeTo(2);
2906 errorsY[i].ResizeTo(2);
2907 covY[i].ResizeTo(2,2);
2908 paramsZ[i].ResizeTo(2);
2909 errorsZ[i].ResizeTo(2);
2910 covZ[i].ResizeTo(2,2);
2911 fittersY[i]->ClearPoints();
2912 fittersZ[i]->ClearPoints();
2918 for (Int_t irow=0;irow<159;irow++) {
2919 AliTPCclusterMI *c=track->GetClusterPointer(irow);
2921 if ((c->GetDetector()%36)!=isec) continue;
2922 Double_t dx = c->GetX()-fXmiddle;
2923 Double_t x[2]={dx, dx*dx};
2924 Double_t yfit = pyf[0]+pyf[1]*dx;
2925 Double_t zfit = pzf[0]+pzf[1]*dx;
2926 Double_t yfitC = pyfc[0]+pyfc[1]*dx;
2927 Double_t zfitC = pzfc[0]+pzfc[1]*dx;
2928 Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
2929 if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
2930 if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
2931 if (dedge<kedgey) continue;
2933 corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
2934 c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
2935 if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
2937 vecX[countRes]=c->GetX();
2938 vecY[countRes]=c->GetY()-yfit;
2939 vecZ[countRes]=c->GetZ()-zfit;
2942 // Fill THnSparse cluster residuals
2943 // use only primary candidates with ITS signal
2944 if (fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){
2945 Double_t resVector[5];
2946 resVector[1]= 9.*gphi/TMath::Pi();
2947 resVector[2]= TMath::Sqrt(c->GetX()*c->GetX()+c->GetY()*c->GetY());
2948 resVector[3]= c->GetZ()/resVector[2];
2951 resVector[0]= c->GetY()-yfitC;
2952 fClusterDelta[0]->Fill(resVector);
2953 resVector[0]= c->GetZ()-zfitC;
2954 fClusterDelta[1]->Fill(resVector);
2956 if (c->GetRow()<kdrow0Fit) continue;
2957 if (c->GetRow()>159-kdrow1Fit) continue;
2960 if (c->GetDetector()>35){
2961 if (c->GetX()<fXquadrant){
2962 if (yfit<-kPRFWidth) fittersY[1]->AddPoint(x,c->GetY(),0.1);
2963 if (yfit<-kPRFWidth) fittersZ[1]->AddPoint(x,c->GetZ(),0.1);
2964 if (yfit>kPRFWidth) fittersY[2]->AddPoint(x,c->GetY(),0.1);
2965 if (yfit>kPRFWidth) fittersZ[2]->AddPoint(x,c->GetZ(),0.1);
2967 if (c->GetX()>fXquadrant){
2968 if (yfit<-kPRFWidth) fittersY[3]->AddPoint(x,c->GetY(),0.1);
2969 if (yfit<-kPRFWidth) fittersZ[3]->AddPoint(x,c->GetZ(),0.1);
2970 if (yfit>kPRFWidth) fittersY[4]->AddPoint(x,c->GetY(),0.1);
2971 if (yfit>kPRFWidth) fittersZ[4]->AddPoint(x,c->GetZ(),0.1);
2974 if (c->GetDetector()<36){
2975 fittersY[0]->AddPoint(x,c->GetY(),0.1);
2976 fittersZ[0]->AddPoint(x,c->GetZ(),0.1);
2982 for (Int_t i=0;i<5;i++) {
2983 npoints[i] = fittersY[i]->GetNpoints();
2984 if (npoints[i]>=kMinClusterQ){
2986 fittersY[i]->Eval();
2987 Double_t chi2FacY = TMath::Sqrt(fittersY[i]->GetChisquare()/(fittersY[i]->GetNpoints()-2));
2989 fittersY[i]->GetParameters(paramsY[i]);
2990 fittersY[i]->GetErrors(errorsY[i]);
2991 fittersY[i]->GetCovarianceMatrix(covY[i]);
2992 // renormalize errors
2993 errorsY[i][0]*=chi2FacY;
2994 errorsY[i][1]*=chi2FacY;
2995 covY[i](0,0)*=chi2FacY*chi2FacY;
2996 covY[i](0,1)*=chi2FacY*chi2FacY;
2997 covY[i](1,0)*=chi2FacY*chi2FacY;
2998 covY[i](1,1)*=chi2FacY*chi2FacY;
3000 fittersZ[i]->Eval();
3001 Double_t chi2FacZ = TMath::Sqrt(fittersZ[i]->GetChisquare()/(fittersZ[i]->GetNpoints()-2));
3003 fittersZ[i]->GetParameters(paramsZ[i]);
3004 fittersZ[i]->GetErrors(errorsZ[i]);
3005 fittersZ[i]->GetCovarianceMatrix(covZ[i]);
3006 // renormalize errors
3007 errorsZ[i][0]*=chi2FacZ;
3008 errorsZ[i][1]*=chi2FacZ;
3009 covZ[i](0,0)*=chi2FacZ*chi2FacZ;
3010 covZ[i](0,1)*=chi2FacZ*chi2FacZ;
3011 covZ[i](1,0)*=chi2FacZ*chi2FacZ;
3012 covZ[i](1,1)*=chi2FacZ*chi2FacZ;
3016 // void UpdateSectorKalman
3018 for (Int_t i0=0;i0<5;i0++){
3019 for (Int_t i1=i0+1;i1<5;i1++){
3020 if(npoints[i0]<kMinClusterQ) continue;
3021 if(npoints[i1]<kMinClusterQ) continue;
3022 TMatrixD v0(4,1),v1(4,1); // measurement
3023 TMatrixD cov0(4,4),cov1(4,4); // covariance
3025 v0(0,0)= paramsY[i0][0]; v1(0,0)= paramsY[i1][0];
3026 v0(1,0)= paramsY[i0][1]; v1(1,0)= paramsY[i1][1];
3027 v0(2,0)= paramsZ[i0][0]; v1(2,0)= paramsZ[i1][0];
3028 v0(3,0)= paramsZ[i0][1]; v1(3,0)= paramsZ[i1][1];
3030 cov0(0,0) = covY[i0](0,0);
3031 cov0(1,1) = covY[i0](1,1);
3032 cov0(1,0) = covY[i0](1,0);
3033 cov0(0,1) = covY[i0](0,1);
3034 cov0(2,2) = covZ[i0](0,0);
3035 cov0(3,3) = covZ[i0](1,1);
3036 cov0(3,2) = covZ[i0](1,0);
3037 cov0(2,3) = covZ[i0](0,1);
3039 cov1(0,0) = covY[i1](0,0);
3040 cov1(1,1) = covY[i1](1,1);
3041 cov1(1,0) = covY[i1](1,0);
3042 cov1(0,1) = covY[i1](0,1);
3043 cov1(2,2) = covZ[i1](0,0);
3044 cov1(3,3) = covZ[i1](1,1);
3045 cov1(3,2) = covZ[i1](1,0);
3046 cov1(2,3) = covZ[i1](0,1);
3050 if (TMath::Abs(pyf[1])<0.8){ //angular cut
3051 UpdateSectorKalman(isec, i0,i1, &v0,&cov0,&v1,&cov1);
3057 // Dump debug information
3059 if (fStreamLevel>0){
3060 // get vertex position
3062 TTreeSRedirector *cstream = GetDebugStreamer();
3066 for (Int_t i0=0;i0<5;i0++){
3067 for (Int_t i1=i0;i1<5;i1++){
3068 if (i0==i1) continue;
3069 if(npoints[i0]<kMinClusterQ) continue;
3070 if(npoints[i1]<kMinClusterQ) continue;
3071 (*cstream)<<"sectorAlign"<<
3072 "run="<<fRun<< // run number
3073 "event="<<fEvent<< // event number
3074 "time="<<fTime<< // time stamp of event
3075 "trigger="<<fTrigger<< // trigger
3076 "triggerClass="<<&fTriggerClass<< // trigger
3077 "mag="<<fMagF<< // magnetic field
3078 "isec="<<isec<< // current sector
3080 "xref="<<fXmiddle<< // reference X
3081 "vPosG.="<<&vPosG<< // vertex position in global system
3082 "vPosL.="<<&vPosL<< // vertex position in local system
3083 "vImpact.="<<&vImpact<< // track impact parameter
3084 //"tofSignal="<<tofSignal<< // tof signal
3085 "tpcPosG.="<<&tpcPosG<< // global position of track at the middle of fXmiddle
3086 "lphi="<<lphi<< // expected local phi angle - if vertex at 0
3087 "gphi="<<gphi<< // expected global phi if vertex at 0
3089 "kyE="<<kyE<< // expect ky - assiming pirmary track
3090 "kzE="<<kzE<< // expected kz - assuming primary tracks
3091 "salpha="<<alpha<< // sector alpha
3092 "scos="<<scos<< // tracking cosinus
3093 "ssin="<<ssin<< // tracking sinus
3097 "nf="<<nf<< // total number of points
3098 "pyf.="<<&pyf<< // full OROC fit y
3099 "pzf.="<<&pzf<< // full OROC fit z
3100 "vX.="<<&vecX<< // x cluster
3101 "vY.="<<&vecY<< // y residual cluster
3102 "vZ.="<<&vecZ<< // z residual cluster
3103 // quadrant and IROC fit
3104 "i0="<<i0<< // quadrant number
3106 "n0="<<npoints[i0]<< // number of points
3107 "n1="<<npoints[i1]<<
3109 "py0.="<<¶msY[i0]<< // parameters
3110 "py1.="<<¶msY[i1]<<
3111 "ey0.="<<&errorsY[i0]<< // errors
3112 "ey1.="<<&errorsY[i1]<<
3113 "chiy0="<<chi2CY[i0]<< // chi2s
3114 "chiy1="<<chi2CY[i1]<<
3116 "pz0.="<<¶msZ[i0]<< // parameters
3117 "pz1.="<<¶msZ[i1]<<
3118 "ez0.="<<&errorsZ[i0]<< // errors
3119 "ez1.="<<&errorsZ[i1]<<
3120 "chiz0="<<chi2CZ[i0]<< // chi2s
3121 "chiz1="<<chi2CZ[i1]<<
3129 void AliTPCcalibAlign::UpdateSectorKalman(Int_t sector, Int_t quadrant0, Int_t quadrant1, TMatrixD *const p0, TMatrixD *const c0, TMatrixD *const p1, TMatrixD *const c1 ){
3132 // tracks are refitted at sector middle
3134 if (fArraySectorIntParam.At(0)==NULL) MakeSectorKalman();
3137 static TMatrixD matHk(4,30); // vector to mesurement
3138 static TMatrixD measR(4,4); // measurement error
3139 // static TMatrixD matQk(2,2); // prediction noise vector
3140 static TMatrixD vecZk(4,1); // measurement
3142 static TMatrixD vecYk(4,1); // Innovation or measurement residual
3143 static TMatrixD matHkT(30,4); // helper matrix Hk transpose
3144 static TMatrixD matSk(4,4); // Innovation (or residual) covariance
3145 static TMatrixD matKk(30,4); // Optimal Kalman gain
3146 static TMatrixD mat1(30,30); // update covariance matrix
3147 static TMatrixD covXk2(30,30); // helper matrix
3149 TMatrixD *vOrig = (TMatrixD*)(fArraySectorIntParam.At(sector));
3150 TMatrixD *cOrig = (TMatrixD*)(fArraySectorIntCovar.At(sector));
3152 TMatrixD vecXk(*vOrig); // X vector
3153 TMatrixD covXk(*cOrig); // X covariance
3157 for (Int_t i=0;i<30;i++)
3158 for (Int_t j=0;j<30;j++){
3160 if (i==j) mat1(i,j)=1;
3164 // matHk - vector to measurement
3166 for (Int_t i=0;i<4;i++)
3167 for (Int_t j=0;j<30;j++){
3177 matHk(0,6*quadrant1+4) = 1.; // delta y
3178 matHk(1,6*quadrant1+0) = 1.; // delta ky
3179 matHk(2,6*quadrant1+5) = 1.; // delta z
3180 matHk(3,6*quadrant1+1) = 1.; // delta kz
3181 // bug fix 24.02 - aware of sign in dx
3182 matHk(0,6*quadrant1+3) = -(*p0)(1,0); // delta x to delta y - through ky
3183 matHk(2,6*quadrant1+3) = -(*p0)(3,0); // delta x to delta z - thorugh kz
3184 matHk(2,6*quadrant1+2) = ((*p0)(0,0)); // y to delta z - through psiz
3186 matHk(0,6*quadrant0+4) = -1.; // delta y
3187 matHk(1,6*quadrant0+0) = -1.; // delta ky
3188 matHk(2,6*quadrant0+5) = -1.; // delta z
3189 matHk(3,6*quadrant0+1) = -1.; // delta kz
3190 // bug fix 24.02 be aware of sign in dx
3191 matHk(0,6*quadrant0+3) = ((*p0)(1,0)); // delta x to delta y - through ky
3192 matHk(2,6*quadrant0+3) = ((*p0)(3,0)); // delta x to delta z - thorugh kz
3193 matHk(2,6*quadrant0+2) = -((*p0)(0,0)); // y to delta z - through psiz
3198 vecZk =(*p1)-(*p0); // measurement
3199 measR =(*c1)+(*c0); // error of measurement
3200 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
3203 matHkT=matHk.T(); matHk.T();
3204 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
3206 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
3207 vecXk += matKk*vecYk; // updated vector
3208 covXk2= (mat1-(matKk*matHk));
3209 covXk = covXk2*covXk;
3216 void AliTPCcalibAlign::MakeSectorKalman(){
3218 // Make a initial Kalman paramaters for IROC - Quadrants alignment
3220 TMatrixD param(5*6,1);
3221 TMatrixD covar(5*6,5*6);
3223 // Set inital parameters
3225 for (Int_t ip=0;ip<5*6;ip++) param(ip,0)=0; // mean alignment to 0
3227 for (Int_t iq=0;iq<5;iq++){
3228 // Initial uncertinty
3229 covar(iq*6+0,iq*6+0) = 0.002*0.002; // 2 mrad
3230 covar(iq*6+1,iq*6+1) = 0.002*0.002; // 2 mrad rotation
3231 covar(iq*6+2,iq*6+2) = 0.002*0.002; // 2 mrad
3233 covar(iq*6+3,iq*6+3) = 0.02*0.02; // 0.2 mm
3234 covar(iq*6+4,iq*6+4) = 0.02*0.02; // 0.2 mm translation
3235 covar(iq*6+5,iq*6+5) = 0.02*0.02; // 0.2 mm
3238 for (Int_t isec=0;isec<36;isec++){
3239 fArraySectorIntParam.AddAt(param.Clone(),isec);
3240 fArraySectorIntCovar.AddAt(covar.Clone(),isec);
3244 void AliTPCcalibAlign::UpdateSectorKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){
3246 // Update kalman vector para0 with vector par1
3249 static TMatrixD matHk(30,30); // vector to mesurement
3250 static TMatrixD measR(30,30); // measurement error
3251 static TMatrixD vecZk(30,1); // measurement
3253 static TMatrixD vecYk(30,1); // Innovation or measurement residual
3254 static TMatrixD matHkT(30,30); // helper matrix Hk transpose
3255 static TMatrixD matSk(30,30); // Innovation (or residual) covariance
3256 static TMatrixD matKk(30,30); // Optimal Kalman gain
3257 static TMatrixD mat1(30,30); // update covariance matrix
3258 static TMatrixD covXk2(30,30); // helper matrix
3260 TMatrixD vecXk(par0); // X vector
3261 TMatrixD covXk(cov0); // X covariance
3266 for (Int_t i=0;i<30;i++)
3267 for (Int_t j=0;j<30;j++){
3269 if (i==j) mat1(i,j)=1;
3271 matHk = mat1; // unit matrix
3273 vecZk = par1; // measurement
3274 measR = cov1; // error of measurement
3275 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
3277 matHkT=matHk.T(); matHk.T();
3278 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
3280 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
3282 vecXk += matKk*vecYk; // updated vector
3283 covXk2= (mat1-(matKk*matHk));
3284 covXk = covXk2*covXk;
3285 CheckCovariance(covXk);
3286 CheckCovariance(cov1);
3288 par0 = vecXk; // update measurement param
3289 cov0 = covXk; // update measurement covar
3292 Double_t AliTPCcalibAlign::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t /*lz*/){
3294 // Get position correction for given sector
3297 TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
3298 if (!param) return 0;
3301 if (lx<fXquadrant) {
3302 if (ly<0) quadrant=1;
3303 if (ly>0) quadrant=2;
3305 if (lx>fXquadrant) {
3306 if (ly<0) quadrant=3;
3307 if (ly>0) quadrant=4;
3310 Double_t a10 = (*param)(quadrant*6+0,0);
3311 Double_t a20 = (*param)(quadrant*6+1,0);
3312 Double_t a21 = (*param)(quadrant*6+2,0);
3313 Double_t dx = (*param)(quadrant*6+3,0);
3314 Double_t dy = (*param)(quadrant*6+4,0);
3315 Double_t dz = (*param)(quadrant*6+5,0);
3316 Double_t deltaX = lx-fXIO;
3317 if (coord==0) return dx;
3318 if (coord==1) return dy+deltaX*a10;
3319 if (coord==2) return dz+deltaX*a20+ly*a21;
3323 Double_t AliTPCcalibAlign::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
3327 if (!Instance()) return 0;
3328 return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz);
3331 void AliTPCcalibAlign::MakeKalman(){
3333 // Make a initial Kalman paramaters for sector Alignemnt
3335 fSectorParamA = new TMatrixD(6*36+6,1);
3336 fSectorParamC = new TMatrixD(6*36+6,1);
3337 fSectorCovarA = new TMatrixD(6*36+6,6*36+6);
3338 fSectorCovarC = new TMatrixD(6*36+6,6*36+6);
3340 // set starting parameters at 0
3342 for (Int_t isec=0;isec<37;isec++)
3343 for (Int_t ipar=0;ipar<6;ipar++){
3344 (*fSectorParamA)(isec*6+ipar,0) =0;
3345 (*fSectorParamC)(isec*6+ipar,0) =0;
3348 // set starting covariance
3350 for (Int_t isec=0;isec<36;isec++)
3351 for (Int_t ipar=0;ipar<6;ipar++){
3353 (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.002*0.002; // 2 mrad
3354 (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.002*0.002;
3357 (*fSectorCovarA)(isec*6+ipar,isec*6+ipar) =0.02*0.02; // 0.2 mm
3358 (*fSectorCovarC)(isec*6+ipar,isec*6+ipar) =0.02*0.02;
3361 (*fSectorCovarA)(36*6+0,36*6+0) =0.04; // common shift y up-up
3362 (*fSectorCovarA)(36*6+1,36*6+1) =0.04; // common shift y down-down
3363 (*fSectorCovarA)(36*6+2,36*6+2) =0.04; // common shift y up-down
3364 (*fSectorCovarA)(36*6+3,36*6+3) =0.004; // common shift phi up-up
3365 (*fSectorCovarA)(36*6+4,36*6+4) =0.004; // common shift phi down-down
3366 (*fSectorCovarA)(36*6+5,36*6+5) =0.004; // common shift phi up-down
3368 (*fSectorCovarC)(36*6+0,36*6+0) =0.04; // common shift y up-up
3369 (*fSectorCovarC)(36*6+1,36*6+1) =0.04; // common shift y down-down
3370 (*fSectorCovarC)(36*6+2,36*6+2) =0.04; // common shift y up-down
3371 (*fSectorCovarC)(36*6+3,36*6+3) =0.004; // common shift phi up-up
3372 (*fSectorCovarC)(36*6+4,36*6+4) =0.004; // common shift phi down-down
3373 (*fSectorCovarC)(36*6+5,36*6+5) =0.004; // common shift phi up-down
3376 void AliTPCcalibAlign::UpdateKalman(Int_t sector0, Int_t sector1, TMatrixD &p0, TMatrixD &c0, TMatrixD &p1, TMatrixD &c1){
3378 // Update Kalman parameters
3379 // Note numbering from 0..36 0..17 IROC 18..35 OROC
3382 if (fSectorParamA==NULL) MakeKalman();
3383 if (CheckCovariance(c0)>0) return;
3384 if (CheckCovariance(c1)>0) return;
3385 const Int_t nelem = 6*36+6;
3388 static TMatrixD matHk(4,nelem); // vector to mesurement
3389 static TMatrixD measR(4,4); // measurement error
3390 static TMatrixD vecZk(4,1); // measurement
3392 static TMatrixD vecYk(4,1); // Innovation or measurement residual
3393 static TMatrixD matHkT(nelem,4); // helper matrix Hk transpose
3394 static TMatrixD matSk(4,4); // Innovation (or residual) covariance
3395 static TMatrixD matKk(nelem,4); // Optimal Kalman gain
3396 static TMatrixD mat1(nelem,nelem); // update covariance matrix
3397 static TMatrixD covXk2(nelem,nelem); // helper matrix
3399 TMatrixD *vOrig = 0;
3400 TMatrixD *cOrig = 0;
3401 vOrig = (sector0%36>=18) ? fSectorParamA:fSectorParamC;
3402 cOrig = (sector0%36>=18) ? fSectorCovarA:fSectorCovarC;
3404 Int_t sec0= sector0%18;
3405 Int_t sec1= sector1%18;
3406 if (sector0>35) sec0+=18;
3407 if (sector1>35) sec1+=18;
3409 TMatrixD vecXk(*vOrig); // X vector
3410 TMatrixD covXk(*cOrig); // X covariance
3414 for (Int_t i=0;i<nelem;i++)
3415 for (Int_t j=0;j<nelem;j++){
3417 if (i==j) mat1(i,j)=1;
3421 // matHk - vector to measurement
3423 for (Int_t i=0;i<4;i++)
3424 for (Int_t j=0;j<nelem;j++){
3434 matHk(0,6*sec1+4) = 1.; // delta y
3435 matHk(1,6*sec1+0) = 1.; // delta ky
3436 matHk(2,6*sec1+5) = 1.; // delta z
3437 matHk(3,6*sec1+1) = 1.; // delta kz
3438 matHk(0,6*sec1+3) = p0(1,0); // delta x to delta y - through ky
3439 matHk(2,6*sec1+3) = p0(3,0); // delta x to delta z - thorugh kz
3440 matHk(2,6*sec1+2) = p0(0,0); // y to delta z - through psiz
3442 matHk(0,6*sec0+4) = -1.; // delta y
3443 matHk(1,6*sec0+0) = -1.; // delta ky
3444 matHk(2,6*sec0+5) = -1.; // delta z
3445 matHk(3,6*sec0+1) = -1.; // delta kz
3446 matHk(0,6*sec0+3) = -p0(1,0); // delta x to delta y - through ky
3447 matHk(2,6*sec0+3) = -p0(3,0); // delta x to delta z - thorugh kz
3448 matHk(2,6*sec0+2) = -p0(0,0); // y to delta z - through psiz
3450 Int_t dsec = (sector1%18)-(sector0%18);
3451 if (dsec<-2) dsec+=18;
3452 if (TMath::Abs(dsec)==1){
3454 // Left right systematic fit part
3457 if (dsec>0) dir= 1.;
3458 if (dsec<0) dir=-1.;
3459 if (sector0>35&§or1>35){
3460 matHk(0,36*6+0)=dir;
3461 matHk(1,36*6+3+0)=dir;
3463 if (sector0<36&§or1<36){
3464 matHk(0,36*6+1)=dir;
3465 matHk(1,36*6+3+1)=dir;
3467 if (sector0<36&§or1>35){
3468 matHk(0,36*6+2)=dir;
3469 matHk(1,36*6+3+2)=dir;
3471 if (sector0>35&§or1<36){
3472 matHk(0,36*6+2)=-dir;
3473 matHk(1,36*6+3+2)=-dir;
3478 vecZk =(p1)-(p0); // measurement
3479 measR =(c1)+(c0); // error of measurement
3480 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
3483 matHkT=matHk.T(); matHk.T();
3484 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
3486 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
3487 vecXk += matKk*vecYk; // updated vector
3488 covXk2= (mat1-(matKk*matHk));
3489 covXk = covXk2*covXk;
3491 if (CheckCovariance(covXk)>0) return;
3500 void AliTPCcalibAlign::UpdateKalman(TMatrixD &par0, TMatrixD &cov0, TMatrixD &par1, TMatrixD &cov1){
3502 // Update kalman vector para0 with vector par1
3505 Int_t nelem = 6*36+6;
3506 static TMatrixD matHk(nelem,nelem); // vector to mesurement
3507 static TMatrixD measR(nelem,nelem); // measurement error
3508 static TMatrixD vecZk(nelem,1); // measurement
3510 static TMatrixD vecYk(nelem,1); // Innovation or measurement residual
3511 static TMatrixD matHkT(nelem,nelem); // helper matrix Hk transpose
3512 static TMatrixD matSk(nelem,nelem); // Innovation (or residual) covariance
3513 static TMatrixD matKk(nelem,nelem); // Optimal Kalman gain
3514 static TMatrixD mat1(nelem,nelem); // update covariance matrix
3515 static TMatrixD covXk2(nelem,nelem); // helper matrix
3517 TMatrixD vecXk(par0); // X vector
3518 TMatrixD covXk(cov0); // X covariance
3523 for (Int_t i=0;i<nelem;i++)
3524 for (Int_t j=0;j<nelem;j++){
3526 if (i==j) mat1(i,j)=1;
3528 matHk = mat1; // unit matrix
3530 vecZk = par1; // measurement
3531 measR = cov1; // error of measurement
3532 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
3534 matHkT=matHk.T(); matHk.T();
3535 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
3537 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
3539 vecXk += matKk*vecYk; // updated vector
3540 covXk2= (mat1-(matKk*matHk));
3541 covXk = covXk2*covXk;
3543 CheckCovariance(cov0);
3544 CheckCovariance(cov1);
3545 CheckCovariance(covXk);
3547 par0 = vecXk; // update measurement param
3548 cov0 = covXk; // update measurement covar
3554 Int_t AliTPCcalibAlign::CheckCovariance(TMatrixD &covar){
3556 // check the consistency of covariance matrix
3558 Int_t ncols = covar.GetNcols();
3559 Int_t nrows= covar.GetNrows();
3560 const Float_t kEpsilon = 0.0001;
3566 printf("Error 0 - wrong matrix\n");
3570 // 1. Check that the non diagonal elements
3572 for (Int_t i0=0;i0<nrows;i0++)
3573 for (Int_t i1=i0+1;i1<ncols;i1++){
3574 Double_t r0 = covar(i0,i1)/TMath::Sqrt(covar(i0,i0)*covar(i1,i1));
3575 Double_t r1 = covar(i1,i0)/TMath::Sqrt(covar(i0,i0)*covar(i1,i1));
3576 if (TMath::Abs(r0-r1)>kEpsilon){
3577 printf("Error 1 - non symetric matrix %d\t%d\t%f",i0,i1,r1-r0);
3580 if (TMath::Abs(r0)>=1){
3581 printf("Error 2 - Wrong correlation %d\t%d\t%f\n",i0,i1,r0);
3584 if (TMath::Abs(r1)>=1){
3585 printf("Error 3 - Wrong correlation %d\t%d\t%f\n",i0,i1,r1);
3594 void AliTPCcalibAlign::MakeReportDy(TFile *output){
3596 // Draw histogram of dY
3598 Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
3599 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
3601 AliTPCcalibAlign *align = this;
3602 TVectorD vecOOP(36);
3603 TVectorD vecOOM(36);
3604 TVectorD vecOIP(36);
3605 TVectorD vecOIM(36);
3606 TVectorD vecOIS(36);
3607 TVectorD vecSec(36);
3608 TCanvas * cOROCdy = new TCanvas("OROC dy","OROC dy",900,600);
3609 cOROCdy->Divide(6,6);
3610 TCanvas * cIROCdy = new TCanvas("IROC dy","IROC dy",900,600);
3611 cIROCdy->Divide(6,6);
3612 TCanvas * cDy = new TCanvas("Dy","Dy",600,700);
3614 for (Int_t isec=0;isec<36;isec++){
3615 Bool_t isDraw=kFALSE;
3617 cOROCdy->cd(isec+1);
3618 Int_t secPlus = (isec%18==17)? isec-17:isec+1;
3619 Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
3620 printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
3621 TH1 * hisOOP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus+36);
3622 TH1 * hisOOM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus+36);
3623 TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec);
3626 hisOIS = (TH1F*)(hisOIS->Clone());
3627 hisOIS->SetDirectory(0);
3628 hisOIS->Scale(1./(hisOIS->GetMaximum()+1));
3629 hisOIS->SetLineColor(kmicolors[0]);
3632 vecOIS(isec)=10*hisOIS->GetMean();
3635 hisOOP = (TH1F*)(hisOOP->Clone());
3636 hisOOP->SetDirectory(0);
3637 hisOOP->Scale(1./(hisOOP->GetMaximum()+1));
3638 hisOOP->SetLineColor(kmicolors[1]);
3639 if (isDraw) hisOOP->Draw("same");
3640 if (!isDraw) {hisOOP->Draw(""); isDraw=kTRUE;}
3641 vecOOP(isec)=10*hisOOP->GetMean();
3644 hisOOM = (TH1F*)(hisOOM->Clone());
3645 hisOOM->SetDirectory(0);
3646 hisOOM->Scale(1/(hisOOM->GetMaximum()+1));
3647 hisOOM->SetLineColor(kmicolors[3]);
3648 if (isDraw) hisOOM->Draw("same");
3649 if (!isDraw) {hisOOM->Draw(""); isDraw=kTRUE;}
3650 vecOOM(isec)=10*hisOOM->GetMean();
3655 for (Int_t isec=0;isec<36;isec++){
3656 Bool_t isDraw=kFALSE;
3657 cIROCdy->cd(isec+1);
3658 Int_t secPlus = (isec%18==17)? isec-17:isec+1;
3659 Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
3660 printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
3661 TH1 * hisOIP= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secPlus);
3662 TH1 * hisOIM= align->GetHisto(AliTPCcalibAlign::kY,isec+36,secMinus);
3663 TH1 * hisOIS= align->GetHisto(AliTPCcalibAlign::kY,isec+36,isec);
3665 hisOIS = (TH1F*)(hisOIS->Clone());
3666 hisOIS->SetDirectory(0);
3667 hisOIS->Scale(1./(hisOIS->GetMaximum()+1));
3668 hisOIS->SetLineColor(kmicolors[0]);
3671 vecOIS(isec)=10*hisOIS->GetMean();
3674 hisOIP = (TH1F*)(hisOIP->Clone());
3675 hisOIP->SetDirectory(0);
3676 hisOIP->Scale(1./(hisOIP->GetMaximum()+1));
3677 hisOIP->SetLineColor(kmicolors[1]);
3678 if (isDraw) hisOIP->Draw("same");
3679 if (!isDraw) {hisOIP->Draw(""); isDraw=kTRUE;}
3680 hisOIP->Draw("same");
3681 vecOIP(isec)=10*hisOIP->GetMean();
3684 hisOIM = (TH1F*)(hisOIM->Clone());
3685 hisOIM->SetDirectory(0);
3686 hisOIM->Scale(1/(hisOIM->GetMaximum()+1));
3687 hisOIM->SetLineColor(kmicolors[3]);
3688 if (isDraw) hisOIM->Draw("same");
3689 if (!isDraw) {hisOIM->Draw(""); isDraw=kTRUE;}
3690 vecOIM(isec)=10*hisOIM->GetMean();
3693 TGraph* grOIM = new TGraph(36,vecSec.GetMatrixArray(),vecOIM.GetMatrixArray());
3694 TGraph* grOIP = new TGraph(36,vecSec.GetMatrixArray(),vecOIP.GetMatrixArray());
3695 TGraph* grOIS = new TGraph(36,vecSec.GetMatrixArray(),vecOIS.GetMatrixArray());
3696 TGraph* grOOM = new TGraph(36,vecSec.GetMatrixArray(),vecOOM.GetMatrixArray());
3697 TGraph* grOOP = new TGraph(36,vecSec.GetMatrixArray(),vecOOP.GetMatrixArray());
3699 grOIS->SetMarkerStyle(kmimarkers[0]);
3700 grOIP->SetMarkerStyle(kmimarkers[1]);
3701 grOIM->SetMarkerStyle(kmimarkers[3]);
3702 grOOP->SetMarkerStyle(kmimarkers[1]);
3703 grOOM->SetMarkerStyle(kmimarkers[3]);
3704 grOIS->SetMarkerColor(kmicolors[0]);
3705 grOIP->SetMarkerColor(kmicolors[1]);
3706 grOIM->SetMarkerColor(kmicolors[3]);
3707 grOOP->SetMarkerColor(kmicolors[1]);
3708 grOOM->SetMarkerColor(kmicolors[3]);
3709 grOIS->SetLineColor(kmicolors[0]);
3710 grOIP->SetLineColor(kmicolors[1]);
3711 grOIM->SetLineColor(kmicolors[3]);
3712 grOOP->SetLineColor(kmicolors[1]);
3713 grOOM->SetLineColor(kmicolors[3]);
3714 grOIS->SetMaximum(1.5);
3715 grOIS->SetMinimum(-1.5);
3716 grOIS->GetXaxis()->SetTitle("Sector number");
3717 grOIS->GetYaxis()->SetTitle("#Delta_{y} (mm)");
3727 cOROCdy->SaveAs("picAlign/OROCOROCdy.eps");
3728 cOROCdy->SaveAs("picAlign/OROCOROCdy.gif");
3729 cOROCdy->SaveAs("picAlign/OROCOROCdy.root");
3731 cIROCdy->SaveAs("picAlign/OROCIROCdy.eps");
3732 cIROCdy->SaveAs("picAlign/OROCIROCdy.gif");
3733 cIROCdy->SaveAs("picAlign/OROCIROCdy.root");
3735 cDy->SaveAs("picAlign/Sectordy.eps");
3736 cDy->SaveAs("picAlign/Sectordy.gif");
3737 cDy->SaveAs("picAlign/Sectordy.root");
3740 cOROCdy->Write("OROCOROCDy");
3741 cIROCdy->Write("OROCIROCDy");
3742 cDy->Write("SectorDy");
3746 void AliTPCcalibAlign::MakeReportDyPhi(TFile */*output*/){
3750 Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
3751 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
3753 AliTPCcalibAlign *align = this;
3754 TCanvas * cOROCdyPhi = new TCanvas("OROC dyphi","OROC dyphi",900,600);
3755 cOROCdyPhi->Divide(6,6);
3756 for (Int_t isec=0;isec<36;isec++){
3757 cOROCdyPhi->cd(isec+1);
3758 Int_t secPlus = (isec%18==17)? isec-17:isec+1;
3759 Int_t secMinus= (isec%18==0) ? isec+17:isec-1;
3760 //printf("%d\t%d\t%d\n",isec,secPlus,secMinus);
3762 TProfile * profdyphiOOP=0,*profdyphiOOM=0,*profdyphiOOS=0;
3763 htemp = (TH2F*) (align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secPlus+36));
3764 if (htemp) profdyphiOOP= htemp->ProfileX();
3765 htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,secMinus+36));
3766 if (htemp) profdyphiOOM= htemp->ProfileX();
3767 htemp = (TH2F*)(align->GetHisto(AliTPCcalibAlign::kYPhi,isec+36,isec));
3768 if (htemp) profdyphiOOS= htemp->ProfileX();
3771 profdyphiOOS->SetLineColor(kmicolors[0]);
3772 profdyphiOOS->SetMarkerStyle(kmimarkers[0]);
3773 profdyphiOOS->SetMarkerSize(0.2);
3774 profdyphiOOS->SetMaximum(0.5);
3775 profdyphiOOS->SetMinimum(-0.5);
3776 profdyphiOOS->SetXTitle("tan(#phi)");
3777 profdyphiOOS->SetYTitle("#DeltaY (cm)");
3780 profdyphiOOP->SetLineColor(kmicolors[1]);
3781 profdyphiOOP->SetMarkerStyle(kmimarkers[1]);
3782 profdyphiOOP->SetMarkerSize(0.2);
3785 profdyphiOOM->SetLineColor(kmicolors[3]);
3786 profdyphiOOM->SetMarkerStyle(kmimarkers[3]);
3787 profdyphiOOM->SetMarkerSize(0.2);
3790 profdyphiOOS->Draw();
3792 if (profdyphiOOM) profdyphiOOM->Draw("");
3793 if (profdyphiOOP) profdyphiOOP->Draw("");
3795 if (profdyphiOOM) profdyphiOOM->Draw("same");
3796 if (profdyphiOOP) profdyphiOOP->Draw("same");