1 /// \file CalibAlignKalman.C
4 /// gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
5 /// gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
6 /// gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
8 /// gROOT->LoadMacro("$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C+");
10 /// AliTPCTransformation::BuildBasicFormulas();
12 /// AliXRDPROOFtoolkit tool;
13 /// chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0,50000);
14 /// chainPoints->Lookup();
16 /// chainMS = tool.MakeChainRandom("kalmanFit.list","kf",0,50000);
17 /// chainMS->Lookup();
19 /// chainFP = tool.MakeChainRandom("kalmanFit.list","filter",0,50000);
20 /// chainFP->Lookup();
22 /// CalibAlignKalmanFit(40000,1);
23 /// kalmanFit0->DumpCorelation(0.8);
24 /// TFile f("kalmanfitTPC.root");
38 #include "TLinearFitter.h"
42 #include "TEntryList.h"
44 #include "THnSparse.h"
47 #include "AliSysInfo.h"
48 #include "AliExternalTrackParam.h"
49 #include "TTreeStream.h"
50 #include "AliTrackPointArray.h"
52 #include "AliTPCTransformation.h"
53 #include "AliTPCkalmanFit.h"
54 #include "AliMathBase.h"
55 #include "AliXRDPROOFtoolkit.h"
63 const Float_t krmsYcut = 0.2; // cluster RMS cut in Y - residuals form local fit
64 const Float_t krmsZcut = 0.2; // cluster RMS cut in Z - residuals from local fit
65 const Float_t krmsYcutGlobal= 0.2; // cluster RMS cut in Y - residuals form global fit
66 const Float_t krmsZcutGlobal= 2.0; // cluster RMS cut in Z - residuals from global fit
67 const Float_t kSigmaCut = 5.; // clusters to be removed
68 const Int_t knclCut = 80; // minimal number of clusters
69 const Double_t kArmCut = 50.; // minimal level arm
70 const Double_t kNsigma = 5.; // cut on track level
72 // mult. scattering cuts
74 TCut cutClY("rms09.fElements[0]<0.15");
75 TCut cuttY("rms09.fElements[1]/rms09.fElements[0]<0.9");
76 TCut cutkY("rms09.fElements[2]<0.015");
77 TCut cutClZ("rms09.fElements[3]<0.2");
78 TCut cuttZ("rms09.fElements[4]/rms09.fElements[3]<0.9");
79 TCut cutkZ("rms09.fElements[5]<0.015");
80 TCut cutMS=cutClY+cuttY+cutkY+cutClZ+cuttZ+cutkZ
82 TMatrixD cutMatrix(4*7,2);
83 const Double_t rmsCut09[6]={0.15,0.9,0.015, 0.2, 0.9, 0.015};
90 Int_t toSkipOffset = 0;
91 Int_t toSkipTrack = 2;
92 Int_t toSkipTrackOffset = 0;
93 Int_t isFilterTest = 0;
97 TCut * cSide[4]={0,0,0,0};
98 TCut *cP3[4]={0,0,0,0};
99 TCut *cSP3[4]={0,0,0,0};
100 TCut *cP4[4]={0,0,0,0};
101 TCut *cM4[4]={0,0,0,0};
102 TCut *cA[4]={0,0,0,0};
110 TChain *chainPoints=0;
112 AliTPCkalmanFit * kalmanFitNew=0;
113 AliTPCkalmanFit * kalmanFitOrig=0;
114 AliTPCkalmanFit * kalmanFitApply=0;
116 AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump);
118 AliTPCkalmanFit * SetupFit();
120 void AddFitFieldCage(AliTPCkalmanFit *kalmanFit);
121 void AddPhiScaling(AliTPCkalmanFit *kalmanFit);
122 void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin );
123 void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
124 void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit);
125 void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit);
126 void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
127 void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
128 void AddDrift(AliTPCkalmanFit *kalmanFit);
131 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream);
133 TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream);
134 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset);
136 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump);
138 void CalibAlignKalman(Int_t npoints, Int_t maxFiles, Int_t startFile, Int_t trackDump, Int_t nSkipTrack, Int_t nSkipTrackOffset, Int_t nSkip, Int_t nSkipOffset, Int_t bfilterTest){
141 AliTPCTransformation::BuildBasicFormulas();
143 toSkipOffset= nSkipOffset;
144 toSkipTrack = nSkipTrack;
145 toSkipTrackOffset = nSkipTrackOffset;
146 isFilterTest = bfilterTest;
148 // read the transformation to be applied
149 TFile ftrafo("kalmanFitApply.root");
150 kalmanFitApply = (AliTPCkalmanFit *)ftrafo.Get("kalmanFitNew");
151 if (kalmanFitApply) {
152 printf("Loaded transforamtion\n");
153 kalmanFitApply->DumpCalib("IROCOROC");
154 kalmanFitApply->InitTransformation();
156 printf("Not trnasformation specified\n");
161 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
162 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
164 AliXRDPROOFtoolkit tool;
165 chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0, maxFiles, startFile);
166 CalibAlignKalmanFit(npoints, trackDump);
172 AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump){
173 /// Fitting procedure
175 AliTPCTransformation::BuildBasicFormulas();
177 kalmanFitNew = SetupFit();
178 kalmanFitOrig = SetupFit();
179 kalmanFitNew->Init();
180 kalmanFitOrig->Init();
181 return FitPointsLinear(maxTracks,trackDump);
186 AliTPCkalmanFit * SetupFit(){
189 AliTPCkalmanFit *kalmanFit = new AliTPCkalmanFit;
190 AddFitFieldCage(kalmanFit);
191 AddPhiScaling(kalmanFit);
193 AddZShift(kalmanFit,3,3);
194 AddZTilting(kalmanFit,3,3);
195 // AddLocalXYMisalignment(kalmanFit);
196 // AddLocalXYMisalignmentSector(kalmanFit);
197 AddAlignSectorFourier(kalmanFit,4,4);
198 AddAlignOROCIROCFourier(kalmanFit,5,5);
207 cSide[0] = new TCut("cutAA","p0In.fP[1]>0&&p1In.fP[1]>0");
208 cSide[1] = new TCut("cutCC","p0In.fP[1]<0&&p1In.fP[1]<0");
209 cSide[2] = new TCut("cutAC","p0In.fP[1]>0&&p1In.fP[1]<0");
210 cSide[3] = new TCut("cutCA","p0In.fP[1]<0&&p1In.fP[1]>0");
212 TH1F * phisP3 = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
213 TH1F * phisSP3 = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
214 TH1F * phisP4 = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
215 TH1F * phisM4 = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
219 TF1 *fg = new TF1("fg","gaus");
220 for (Int_t iter=0; iter<2; iter++){
221 for (Int_t ic=0;ic<4;ic++){
229 chainPoints->Draw("p0In.fP[3]+p1In.fP[3]>>hhisP3",*cA[ic],"goff");
230 phisP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisP3->GetMaximumBin())-0.003,phisP3->GetBinCenter(phisP3->GetMaximumBin())+0.003);
231 cP3[ic]=new TCut(Form("abs(p0In.fP[3]+p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
232 cutMatrix(7*ic+0,0) = fg->GetParameter(1);
233 cutMatrix(7*ic+0,1) = fg->GetParameter(2);
237 chainPoints->Draw("p0Out.fP[3]-p0In.fP[3]>>hhisSP3",*cA[ic],"goff");
238 phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
239 cSP3[ic]=new TCut(Form("abs(p0Out.fP[3]-p0In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
240 cutMatrix(7*ic+1,0) = fg->GetParameter(1);
241 cutMatrix(7*ic+1,1) = fg->GetParameter(2);
243 chainPoints->Draw("p1Out.fP[3]-p1In.fP[3]>>hhisSP3",*cA[ic],"goff");
244 phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
245 *cSP3[ic]+=Form("abs(p1Out.fP[3]-p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
246 cutMatrix(7*ic+2,0) = fg->GetParameter(1);
247 cutMatrix(7*ic+2,1) = fg->GetParameter(2);
251 chainPoints->Draw("p0Out.fP[4]>>hhisP4",*cA[ic],"goff");
252 phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
253 cP4[ic]=new TCut(Form("abs(p0Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
254 cutMatrix(7*ic+3,0) = fg->GetParameter(1);
255 cutMatrix(7*ic+3,1) = fg->GetParameter(2);
256 chainPoints->Draw("p1Out.fP[4]>>hhisP4",*cA[ic],"goff");
257 phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
258 *cP4[ic]+=Form("abs(p1Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
259 cutMatrix(7*ic+4,0) = fg->GetParameter(1);
260 cutMatrix(7*ic+4,1) = fg->GetParameter(2);
265 chainPoints->Draw("p0Out.fP[4]-p0In.fP[4]>>hhisM4",*cA[ic],"goff");
266 phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
267 cM4[ic]=new TCut(Form("abs(p0Out.fP[4]-p0In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
268 cutMatrix(7*ic+5,0) = fg->GetParameter(1);
269 cutMatrix(7*ic+5,1) = fg->GetParameter(2);
271 chainPoints->Draw("p1Out.fP[4]-p1In.fP[4]>>hhisM4",*cA[ic],"goff");
272 phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
273 *cM4[ic]+=Form("abs(p1Out.fP[4]-p1In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
274 cutMatrix(7*ic+6,0) = fg->GetParameter(1);
275 cutMatrix(7*ic+6,1) = fg->GetParameter(2);
280 *cA[ic]= *cSide[ic]+*cP3[ic]+*cSP3[ic]+*cP4[ic]+*cM4[ic];
285 cutAll = (*cA[0])||(*cA[1])||(*cA[2])||(*cA[3])+"abs(mag)<0.01&&ncont>0&&p.fNPoints>120";
287 delete phisP3; // = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
288 delete phisSP3; // = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
289 delete phisP4;// = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
290 delete phisM4;// = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
294 chainPoints->Draw(">>listEL",cutAll,"entryList");
295 elist = (TEntryList*)gDirectory->Get("listEL");
296 chainPoints->SetEntryList(elist);
297 elist->SetDirectory(0);
302 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump){
303 // create debug streamers
305 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPC.root");
306 TTreeSRedirector *pcstreamOrig = new TTreeSRedirector("kalmanfitTPCOrig.root");
307 pcstream->GetFile()->cd();
308 cutMatrix.Write("cutMarix");
309 elist->Write("eventList");
312 AliTrackPointArray *pointsNS=0;
315 AliExternalTrackParam *param0=0;
316 AliExternalTrackParam *param1=0;
317 chainPoints->SetBranchAddress("p.",&pointsNS);
318 chainPoints->SetBranchAddress("p0In.",¶m0);
319 chainPoints->SetBranchAddress("p1In.",¶m1);
320 chainPoints->SetBranchAddress("mag",&mag);
321 chainPoints->SetBranchAddress("time",&time);
323 printf("\n*\n*\n*Selected entries = %d\n*\n*\n*",Int_t(elist->GetN()));
326 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
327 if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
328 Int_t entry=chainPoints->GetEntryNumber(itrack);
329 chainPoints->GetEntry(entry);
330 if (accepted>maxTracks) break;
332 AliTrackPointArray *points = AliTPCkalmanFit::SortPoints(*pointsNS);
333 if (kalmanFitApply) kalmanFitApply->ApplyCalibration(points,-1.);
335 // estimate and filter scattering
337 TVectorD *vecRMS09 = EstimateScatteringKalmanLinear(*points,*param0,*param1,pcstream);
338 if (!vecRMS09) continue;
340 if ((*vecRMS09)[0] >rmsCut09[0]) isOK=kFALSE;
341 if ((*vecRMS09)[1]/(*vecRMS09)[0] >rmsCut09[1]) isOK=kFALSE;;
342 if ((*vecRMS09)[2] >rmsCut09[2]) isOK=kFALSE;
343 if ((*vecRMS09)[3] >rmsCut09[3]) isOK=kFALSE;
344 if ((*vecRMS09)[4]/(*vecRMS09)[0] >rmsCut09[4]) isOK=kFALSE;
345 if ((*vecRMS09)[5] >rmsCut09[5]) isOK=kFALSE;
346 if (!isOK || isFilterTest) {
350 kalmanFitNew->PropagateTime(time);
353 for (Int_t idir=-1; idir<=1; idir++){
354 AliTrackPointArray *pointsF = FilterPoints(*points,idir, pcstream);
355 if (!pointsF) continue;
356 AliTrackPointArray *spointsF = 0;
357 // we skip points for alignemnt but not for QA
358 if (idir==0) spointsF = SkipPoints(*pointsF, toSkip*2, toSkipOffset);
359 if (idir!=0) spointsF = SkipPoints(*pointsF, toSkip, toSkipOffset);
360 if (idir==0) accepted++;
362 if (accepted%50==0) {
363 kalmanFitNew->FitTrackLinear(*pointsF, pcstream);
365 if (idir==0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
366 if (idir!=0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
368 if (idir==0) kalmanFitNew->DumpTrackLinear(*pointsF,pcstream);
369 if (idir==0) kalmanFitOrig->DumpTrackLinear(*pointsF,pcstreamOrig);
370 if (accepted%trackDump==0) {
371 printf("%d\n", accepted);
373 AliSysInfo::AddStamp("trackFit", accepted,itrack);
379 pcstream->GetFile()->cd();
380 kalmanFitNew->Write("kalmanFit");
381 pcstreamOrig->GetFile()->cd();
382 kalmanFitOrig->Write("kalmanFitOrig");
383 pcstreamOrig->GetFile()->cd();
384 if (kalmanFitApply) kalmanFitApply->Write("kalmanFitApply");
391 void QAPointsLinear(Int_t maxTracks, Int_t trackDump){
392 /// check the consistency of kalman fit
393 /// Apply transformation
395 // create debug streeamers
396 TTreeSRedirector *pcstreamNonCalib = new TTreeSRedirector("kalmanfitTPCQANonCalib.root");
397 TTreeSRedirector *pcstreamCalib = new TTreeSRedirector("kalmanfitTPCQACalib.root");
398 TTreeSRedirector *pcstream=0;
399 AliTPCkalmanFit *kalmanFitters[6]={0,0,0,0,0,0};
400 for (Int_t i=0;i<6;i++){
401 kalmanFitters[i]=SetupFit();
404 AliTrackPointArray *points=0;
405 AliExternalTrackParam *param0=0;
406 AliExternalTrackParam *param1=0;
409 chainPoints->SetBranchAddress("p.",&points);
410 chainPoints->SetBranchAddress("mag",&mag);
411 chainPoints->SetBranchAddress("time",&time);
412 chainPoints->SetBranchAddress("p0In.",¶m0);
413 chainPoints->SetBranchAddress("p1In.",¶m1);
419 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
420 if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
421 Int_t entry=chainPoints->GetEntryNumber(itrack);
422 chainPoints->GetEntry(entry);
423 if (accepted>maxTracks) break;
426 AliTrackPointArray pointsCalib(*points);
427 for (Int_t iscalib=0; iscalib<1;iscalib++){
428 if (iscalib>0) kalmanFitNew->ApplyCalibration(&pointsCalib,-1.);
429 if (iscalib==0) pcstream=pcstreamNonCalib;
430 if (iscalib>0) pcstream=pcstreamCalib;
431 for (Int_t idir=-1; idir<=1; idir++){
432 AliTrackPointArray *pointsF = FilterPoints(pointsCalib,idir, pcstream);
433 if (!pointsF) continue;
435 if (idir==0) accepted++;
436 kalmanFitters[iscalib*3+idir+1]->DumpTrackLinear(*pointsF,0);
437 EstimateScatteringKalmanLinear(*pointsF,*param0,*param1,pcstream);
441 if (accepted%trackDump==0) {
442 printf("%d\n", accepted);
445 pcstreamCalib->GetFile()->cd();
446 kalmanFitters[0]->Write("fitUpNonCalib");
447 kalmanFitters[1]->Write("fitUpDownNonCalib");
448 kalmanFitters[2]->Write("fitDownNonCalib");
449 kalmanFitters[3]->Write("fitUpCalib");
450 kalmanFitters[4]->Write("fitUpDownCalib");
451 kalmanFitters[5]->Write("fitDownCalib");
452 delete pcstreamCalib;
453 delete pcstreamNonCalib;
457 void TestScattering(Int_t maxTracks, Int_t trackDump){
458 /// test Multiple scattering algorithm
459 /// Apply transformation
461 /// create debug streeamers
463 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPCMS.root");
466 AliTrackPointArray *points=0;
467 AliExternalTrackParam *param0=0;
468 AliExternalTrackParam *param1=0;
471 chainPoints->SetBranchAddress("p.",&points);
472 chainPoints->SetBranchAddress("mag",&mag);
473 chainPoints->SetBranchAddress("time",&time);
474 chainPoints->SetBranchAddress("p0In.",¶m0);
475 chainPoints->SetBranchAddress("p1In.",¶m1);
478 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
479 Int_t entry=chainPoints->GetEntryNumber(itrack);
480 chainPoints->GetEntry(entry);
481 if (accepted>maxTracks) break;
483 AliTrackPointArray *pointsSorted = AliTPCkalmanFit::SortPoints(*points);
484 EstimateScatteringKalmanLinear(*pointsSorted,*param0,*param1,pcstream);
486 if (accepted%trackDump==0) {
487 printf("%d\n", accepted);
498 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset){
499 /// create new array with skipped points
501 Int_t npoints = points.GetNPoints();
502 Int_t npointsF = (npoints-nskipOffset-1)/nskip;
504 AliTrackPointArray *pointsF= new AliTrackPointArray(npointsF);
506 for (Int_t ipoint=nskipOffset; ipoint<npoints; ipoint+=nskip){
508 if (!points.GetPoint(point,ipoint)) continue;
509 pointsF->AddPoint(used,&point);
511 if (used==npointsF) break;
518 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream){
519 /// Filter points - input points for KalmanFilter
521 TLinearFitter lfitY(2,"pol1");
522 TLinearFitter lfitZ(2,"pol1");
526 lfitY.StoreData(kTRUE);
527 lfitZ.StoreData(kTRUE);
528 Int_t npoints = points.GetNPoints();
529 if (npoints<2) return 0;
530 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
531 Double_t ca = TMath::Cos(currentAlpha);
532 Double_t sa = TMath::Sin(currentAlpha);
534 // 1.b Fit the track in the rotated frame - MakeSeed
536 Double_t maxX =-10000, minX=10000;
537 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
538 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
539 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
540 Double_t rz = points.GetZ()[ipoint];
541 if (dir== 1 && rx<0) continue;
542 if (dir==-1 && rx>0) continue;
543 if (maxX<rx) maxX=rx;
544 if (minX>rx) minX=rx;
545 lfitY.AddPoint(&rx,ry,1);
546 lfitZ.AddPoint(&rx,rz,1);
548 if (TMath::Abs(maxX-minX)<kArmCut) return 0;
549 if (lfitY.GetNpoints()<knclCut) return 0;
553 lfitY.GetParameters(vecY);
554 lfitZ.GetParameters(vecZ);
556 Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
557 Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
558 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
559 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
565 AliTrackPointArray *pointsF=0;
566 for (Int_t iter=0; iter<2;iter++){
567 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
569 if (!points.GetPoint(point,ipoint)) continue;
570 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
571 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
572 Double_t rz = points.GetZ()[ipoint];
573 if (dir== 1 && rx<0) continue;
574 if (dir==-1 && rx>0) continue;
575 Double_t erry = TMath::Sqrt(chi2Y);
576 Double_t errz = TMath::Sqrt(chi2Z);
577 Double_t fy = vecY[0]+vecY[1]*rx;
578 Double_t fz = vecZ[0]+vecZ[1]*rx;
579 if (TMath::Abs(fy-ry)>erry*kSigmaCut) continue;
580 if (TMath::Abs(fz-rz)>errz*kSigmaCut) continue;
582 if (pointsF) pointsF->AddPoint(toBeUsed,&point);
586 (*pcstream)<<"filter"<<
588 "accepted="<<accepted<<
597 if (accepted<knclCut) break;
598 if (iter==0) pointsF = new AliTrackPointArray(toBeUsed);
606 AliTrackPointArray * SortPoints(AliTrackPointArray &points){
607 /// Creates the array - points sorted according radius - neccessay for kalman fit
609 /// 0. choose the frame - rotation angle
611 Int_t npoints = points.GetNPoints();
612 if (npoints<1) return 0;
613 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
614 Double_t ca = TMath::Cos(currentAlpha);
615 Double_t sa = TMath::Sin(currentAlpha);
617 // 1. sort the points
619 Double_t *rxvector = new Double_t[npoints];
620 Int_t *indexes = new Int_t[npoints];
621 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
622 rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
624 TMath::Sort(npoints, rxvector,indexes,kFALSE);
626 AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
627 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
628 if (!points.GetPoint(point,indexes[ipoint])) continue;
629 pointsSorted->AddPoint(ipoint,&point);
636 TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream){
637 /// Algorithm - 0. Fit the track forward and backward
638 /// - 1. Store the current parameters in each point
640 const Int_t kMinPoints= 70;
641 const Double_t kResY = 0.1;
642 const Double_t kResZ = 0.1;
643 const Double_t kMisY = 0.02;
644 const Double_t kMisZ = 0.02;
645 const Double_t kLArm = 120.;
646 const Double_t kACsideFac = 10.;
647 const Double_t kMS0 = kResY/20.;
649 Int_t npoints = points.GetNPoints();
650 if (npoints<kMinPoints) return 0;
652 TVectorD *vecPos[11]={0,0,0,0,0,0,0,0,0,0,0};
653 for (Int_t i=0;i<11;i++){
654 vecPos[i]=new TVectorD(npoints);
658 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
659 Double_t ca = TMath::Cos(currentAlpha);
660 Double_t sa = TMath::Sin(currentAlpha);
663 TMatrixD trParamY(2,1),trCovarY(2,2);
664 TMatrixD trParamZ(2,1),trCovarZ(2,2);
668 TMatrixD matHk(kNmeas,nelem); // vector to mesurement
669 TMatrixD matHkT(nelem,kNmeas); // helper matrix Hk transpose
670 TMatrixD matFk(nelem,nelem);
671 TMatrixD matFkT(nelem,nelem);
672 TMatrixD vecYk(kNmeas,1); // Innovation or measurement residual
673 TMatrixD vecZk(kNmeas,1); // Innovation or measurement residual
674 TMatrixD measR(kNmeas,kNmeas);
675 TMatrixD matSk(kNmeas,kNmeas); // Innovation (or residual) covariance
676 TMatrixD matKk(nelem,kNmeas); // Optimal Kalman gain
677 TMatrixD covXk2(nelem,nelem); // helper matrix
678 TMatrixD covXk3(nelem,nelem); // helper matrix
679 TMatrixD mat1(nelem,nelem);
680 mat1(0,0)=1.; mat1(0,1)=0.;
681 mat1(1,0)=0.; mat1(1,1)=1.;
687 for (Int_t idir=0; idir<2;idir++){
690 for (Int_t ip=0; ip<npoints; ip++){
692 if (idir>0) ipoint = npoints-ip-1;
693 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
694 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
695 Double_t rz = points.GetZ()[ipoint];
696 Int_t volId= points.GetVolumeID()[ipoint];
699 // set initital parameters and covariance - use first and middle point
700 Double_t rxm = ca*points.GetX()[npoints/2]+sa*points.GetY()[npoints/2];
701 Double_t rym = -sa*points.GetX()[npoints/2]+ca*points.GetY()[npoints/2];
702 Double_t rzm = points.GetZ()[npoints/2];
704 trParamY(1,0) = (rym-ry)/(rxm-rx);
706 trParamZ(1,0) = (rzm-rz)/(rxm-rx);
708 trCovarY(0,0) = kResY*kResY;
709 trCovarY(1,1) = (kResY*kResY)/((rxm-rx)*(rxm-rx));
710 trCovarZ(0,0) = kResZ*kResZ;
711 trCovarZ(1,1) = (kResZ*kResZ)/((rxm-rx)*(rxm-rx));
718 if ((volId%36)<18 && (lastVolId%36)>=18){
720 trCovarY(0,0)+=kMisY*kMisY*kACsideFac;
721 trCovarZ(0,0)+=kMisZ*kMisZ*kACsideFac;
722 trCovarY(1,1)+=kACsideFac*(kMisY*kMisY)/(kLArm*kLArm);;
723 trCovarZ(1,1)+=kACsideFac*(kMisZ*kMisZ)/(kLArm*kLArm);
726 if (volId!=lastVolId){
728 trCovarY(0,0)+=kMisY*kMisY;
729 trCovarZ(0,0)+=kMisZ*kMisZ;
730 trCovarY(1,1)+=(kMisY*kMisY)/(kLArm*kLArm);;
731 trCovarZ(1,1)+=(kMisZ*kMisZ)/(kLArm*kLArm);
737 Double_t deltaX=rx-lastX;
738 trParamY(0,0)+=deltaX*trParamY(1,0);
739 trParamZ(0,0)+=deltaX*trParamZ(1,0);
740 matFk(0,0)=1; matFk(0,1)=deltaX;
741 matFk(1,0)=0; matFk(1,1)=1;
742 matFkT=matFk.T(); matFk.T();
743 covXk2=matFk*trCovarY*matFkT;
745 covXk2=matFk*trCovarZ*matFkT;
748 // multiple scattering
749 trCovarY(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
750 trCovarZ(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
755 for (Int_t coord=0; coord<2;coord++){
756 TMatrixD* pvecXk = (coord==0)? &trParamY: &trParamZ;
757 TMatrixD* pcovXk = (coord==0)? &trCovarY: &trCovarZ;
759 TMatrixD& vecXk = *pvecXk;
760 TMatrixD& covXk = *pcovXk;
761 measR(0,0) = (coord==0) ? kResY:kResZ;
762 vecZk(0,0) = (coord==0) ? ry:rz;
764 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
765 matHkT=matHk.T(); matHk.T();
766 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
768 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
770 covXk2= (mat1-(matKk*matHk));
771 covXk3 = covXk2*covXk;
772 vecXk += matKk*vecYk; // updated vector
778 (*vecPos[0])[ipoint]=rx;
779 (*vecPos[1])[ipoint]=ry;
780 (*vecPos[2])[ipoint]=rz;
782 (*vecPos[4*idir+0+3])[ipoint]=trParamY(0,0);
783 (*vecPos[4*idir+1+3])[ipoint]=trParamY(1,0);
785 (*vecPos[4*idir+2+3])[ipoint]=trParamZ(0,0);
786 (*vecPos[4*idir+3+3])[ipoint]=trParamZ(1,0);
792 TVectorD vec(npoints);
794 TVectorD rms09(6); // robust RMS - fraction 0.9
795 TVectorD mean09(6); // robust RMS - fraction 0.9
797 Int_t npoints09 = Int_t(npoints*0.9);
799 vec=(*(vecPos[3])-*(vecPos[1]));
800 rms[0]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms cluster y
801 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
804 vec=(*(vecPos[7])-*(vecPos[3]));
805 rms[1]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track y
806 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
809 vec=(*(vecPos[8])-*(vecPos[4]));
810 rms[2]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track ky
811 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
815 vec=(*(vecPos[5])-*(vecPos[2]));
816 rms[3]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms cluster z
817 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
820 vec=(*(vecPos[9])-*(vecPos[5]));
821 rms[4]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track z
822 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
825 vec=(*(vecPos[10])-*(vecPos[6]));
826 rms[5]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track kz
827 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
838 "mean09.="<<&mean09<<
844 for (Int_t i=0;i<11;i++){
847 return new TVectorD(rms09);
858 void AddFitFieldCage(AliTPCkalmanFit *kalmanFit){
859 /// Add radial scaling due field cage
862 AliTPCTransformation * transformation=0;
865 // linear R scaling and shift
867 for (Int_t iside=0; iside<=1; iside++)
868 for (Int_t ipolR=0; ipolR<2; ipolR++){
869 for (Int_t ipolZ=0; ipolZ<3; ipolZ++){
872 if (ipolR+ipolZ==0) continue;
873 sprintf(tname,"tTPCscalingRPolR%dDr%dSide%d",ipolR,ipolZ,iside);
874 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0, 1);
875 transformation->SetParams(0,0.2,0,&fpar);
876 kalmanFit->AddCalibration(transformation);
883 for (Int_t iside=0; iside<=1; iside++)
884 for (Int_t ipol=0; ipol<3; ipol++){
886 sprintf(tname,"tTPCscalingRIFC%dSide%d",ipol,iside);
887 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRIFC",0,0, 1);
888 transformation->SetParams(0,0.2,0,&fpar);
889 kalmanFit->AddCalibration(transformation);
894 for (Int_t iside=0; iside<=1; iside++)
895 for (Int_t ipol=0; ipol<3; ipol++){
898 sprintf(tname,"tTPCscalingROFC%dSide%d",ipol,iside);
899 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingROFC",0,0, 1);
900 transformation->SetParams(0,0.2,0,&fpar);
901 kalmanFit->AddCalibration(transformation);
906 void AddPhiScaling(AliTPCkalmanFit *kalmanFit){
907 /// Add linear local phi scaling -
908 /// separate IROC/OROC - A side/C side
909 /// "tscalingLocalPhiIROC"
910 /// "tscalingLocalPhiOROC"
914 for (Int_t i=0;i<72;i++){
923 AliTPCTransformation * transformation=0;
925 transformation = new AliTPCTransformation("tscalingLocalPhiIROC", new TBits(maskInner), 0,"TPCscalingPhiLocal",0, 1);
926 transformation->SetParams(0,0.02,0,&fpar);
927 kalmanFit->AddCalibration(transformation);
928 transformation = new AliTPCTransformation("tscalingLocalPhiOROC", new TBits(maskOuter), 0,"TPCscalingPhiLocal",0, 1);
929 transformation->SetParams(0,0.02,0,&fpar);
930 kalmanFit->AddCalibration(transformation);
934 void AddDrift(AliTPCkalmanFit *kalmanFit){
935 /// Add drift velocity transformation
938 AliTPCTransformation * transformation=0;
940 transformation = new AliTPCTransformation("tTPCscalingZDrift", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDr", 0);
941 transformation->SetParams(0,4.,1.,&fpar);
942 kalmanFit->AddCalibration(transformation);
944 transformation = new AliTPCTransformation("tTPCscalingZDriftGy", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDrGy", 0);
945 transformation->SetParams(0,0.2,0.0,&fpar);
946 kalmanFit->AddCalibration(transformation);
952 void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
956 fpar[0]=0; fpar[1]=0; fpar[2]=0;
958 AliTPCTransformation * transformation=0;
962 for (Int_t i=0; i<=ncos;i++){
963 fpar[0]=i; // cos frequency
964 fpar[1]=0; // no sinus
965 fpar[2]=1; // relative misalignment
967 sprintf(tname,"tTPCDeltaZIROCOROCSideA_Cos%d",i);
968 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
969 transformation->SetParams(0,0.03,0,&fpar);
970 kalmanFit->AddCalibration(transformation);
972 sprintf(tname,"tTPCDeltaZIROCOROCSideC_Cos%d",i);
973 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
974 transformation->SetParams(0,0.03,0,&fpar);
975 kalmanFit->AddCalibration(transformation);
979 fpar[2]=0; // absolute misalignment
980 sprintf(tname,"tTPCDeltaZSectorSideA_Cos%d",i);
981 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
982 transformation->SetParams(0,0.1,0,&fpar);
983 kalmanFit->AddCalibration(transformation);
985 sprintf(tname,"tTPCDeltaZSectorSideC_Cos%d",i);
986 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
987 transformation->SetParams(0,0.1,0,&fpar);
988 kalmanFit->AddCalibration(transformation);
991 for (Int_t i=1; i<=nsin;i++){
992 fpar[0]=0; // cos frequency
993 fpar[1]=i; // sinus frequncy
994 fpar[2]=1; // relative misalignment
996 sprintf(tname,"tTPCDeltaZIROCOROCSideA_Sin%d",i);
997 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
998 transformation->SetParams(0,0.03,0,&fpar);
999 kalmanFit->AddCalibration(transformation);
1001 sprintf(tname,"tTPCDeltaZIROCOROCSideC_Sin%d",i);
1002 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1003 transformation->SetParams(0,0.03,0,&fpar);
1004 kalmanFit->AddCalibration(transformation);
1008 fpar[2]=0; // absolute misalignment
1009 sprintf(tname,"tTPCDeltaZSectorSideA_Sin%d",i);
1010 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
1011 transformation->SetParams(0,0.1,0,&fpar);
1012 kalmanFit->AddCalibration(transformation);
1014 sprintf(tname,"tTPCDeltaZSectorSideC_Sin%d",i);
1015 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1016 transformation->SetParams(0,0.1,0,&fpar);
1017 kalmanFit->AddCalibration(transformation);
1025 void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1026 /// z tilting absolute (sector) and relative (IROC-OROC)
1029 fpar[0]=0; fpar[1]=0; fpar[2]=0;
1031 AliTPCTransformation * transformation=0;
1035 for (Int_t i=0; i<=ncos;i++){
1036 fpar[0]=i; // cos frequency
1037 fpar[1]=0; // sinus frequency
1038 fpar[2]=1; // relative misalignment
1040 sprintf(tname,"tTPCTiltingZIROCOROCSideA_Cos%d",i);
1041 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1042 transformation->SetParams(0,0.03,0,&fpar);
1043 kalmanFit->AddCalibration(transformation);
1045 sprintf(tname,"tTPCTiltingZIROCOROCSideC_Cos%d",i);
1046 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1047 transformation->SetParams(0,0.03,0,&fpar);
1048 kalmanFit->AddCalibration(transformation);
1052 fpar[2]=0; // absolute misalignment
1053 sprintf(tname,"tTPCTiltingZSectorSideA_Cos%d",i);
1054 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1055 transformation->SetParams(0,0.1,0,&fpar);
1056 kalmanFit->AddCalibration(transformation);
1058 sprintf(tname,"tTPCTiltingZSectorSideC_Cos%d",i);
1059 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1060 transformation->SetParams(0,0.1,0,&fpar);
1061 kalmanFit->AddCalibration(transformation);
1064 for (Int_t i=1; i<=nsin;i++){
1065 fpar[0]=0; // cos frequency
1066 fpar[1]=i; // sinus frequncy
1067 fpar[2]=1; // relative misalignment
1069 sprintf(tname,"tTPCTiltingZIROCOROCSideA_Sin%d",i);
1070 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1071 transformation->SetParams(0,0.03,0,&fpar);
1072 kalmanFit->AddCalibration(transformation);
1074 sprintf(tname,"tTPCTiltingZIROCOROCSideC_Sin%d",i);
1075 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1076 transformation->SetParams(0,0.03,0,&fpar);
1077 kalmanFit->AddCalibration(transformation);
1081 fpar[2]=0; // absolute misalignment
1082 sprintf(tname,"tTPCTiltingZSectorSideA_Sin%d",i);
1083 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1084 transformation->SetParams(0,0.1,0,&fpar);
1085 kalmanFit->AddCalibration(transformation);
1087 sprintf(tname,"tTPCTiltingZSectorSideC_Sin%d",i);
1088 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1089 transformation->SetParams(0,0.1,0,&fpar);
1090 kalmanFit->AddCalibration(transformation);
1096 void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit){
1100 AliTPCTransformation * transformation=0;
1101 TBits maskInnerA(72);
1102 TBits maskInnerC(72);
1103 for (Int_t i=0;i<72;i++){
1105 if (i%36<18) maskInnerA[i]=kTRUE;
1106 if (i%36>=18) maskInnerC[i]=kTRUE;
1111 transformation = new AliTPCTransformation("tTPCDeltaLxIROCA", new TBits(maskInnerA), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1112 transformation->SetParams(0,0.2,0,&fpar);
1113 kalmanFit->AddCalibration(transformation);
1114 transformation = new AliTPCTransformation("tTPCDeltaLxIROCC", new TBits(maskInnerC), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1115 transformation->SetParams(0,0.2,0,&fpar);
1116 kalmanFit->AddCalibration(transformation);
1118 transformation = new AliTPCTransformation("tTPCDeltaLyIROCA", new TBits(maskInnerA), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1119 transformation->SetParams(0,0.2,0,&fpar);
1120 kalmanFit->AddCalibration(transformation);
1121 transformation = new AliTPCTransformation("tTPCDeltaLyIROCC", new TBits(maskInnerC), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1122 transformation->SetParams(0,0.2,0,&fpar);
1123 kalmanFit->AddCalibration(transformation);
1126 void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit){
1130 AliTPCTransformation * transformation=0;
1133 for (Int_t isec=0; isec<36;isec++){
1138 sprintf(tname,"tTPClocalLxIROC%d",isec);
1139 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1140 transformation->SetParams(0,0.2,0,&fpar);
1141 kalmanFit->AddCalibration(transformation);
1143 sprintf(tname,"tTPClocalLyIROC%d",isec);
1144 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1145 transformation->SetParams(0,0.2,0,&fpar);
1146 kalmanFit->AddCalibration(transformation);
1148 sprintf(tname,"tTPClocalRzIROC%d",isec);
1149 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0, 0);
1150 transformation->SetParams(0,0.2,0,&fpar);
1151 kalmanFit->AddCalibration(transformation);
1155 for (Int_t isec=0; isec<36;isec++){
1156 if (isec%18==fixSector) continue;
1159 mask[isec+36]=kTRUE;
1162 sprintf(tname,"tTPClocalLxSector%d",isec);
1163 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1164 transformation->SetParams(0,0.2,0,&fpar);
1165 kalmanFit->AddCalibration(transformation);
1167 sprintf(tname,"tTPClocalLySector%d",isec);
1168 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1169 transformation->SetParams(0,0.2,0,&fpar);
1170 kalmanFit->AddCalibration(transformation);
1172 sprintf(tname,"tTPClocalRzSector%d",isec);
1173 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0, 0);
1174 transformation->SetParams(0,0.2,0,&fpar);
1175 kalmanFit->AddCalibration(transformation);
1182 void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1186 AliTPCTransformation * transformation=0;
1188 for (Int_t i=0; i<=ncos;i++){
1190 fpar[0]=i; // cos frequency
1191 fpar[1]=0; // no sinus
1192 fpar[2]=1; // relative misalignment
1196 sprintf(tname,"tTPClocalLxIROCOROCSideA_Cos%d",i);
1197 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1198 transformation->SetParams(0,0.03,0,&fpar);
1199 kalmanFit->AddCalibration(transformation);
1200 sprintf(tname,"tTPClocalLxIROCOROCSideC_Cos%d",i);
1201 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1202 transformation->SetParams(0,0.03,0,&fpar);
1203 kalmanFit->AddCalibration(transformation);
1208 sprintf(tname,"tTPClocalLyIROCOROCSideA_Cos%d",i);
1209 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1210 transformation->SetParams(0,0.03,0,&fpar);
1211 kalmanFit->AddCalibration(transformation);
1212 sprintf(tname,"tTPClocalLyIROCOROCSideC_Cos%d",i);
1213 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1214 transformation->SetParams(0,0.03,0,&fpar);
1215 kalmanFit->AddCalibration(transformation);
1221 sprintf(tname,"tTPClocalRzIROCOROCSideA_Cos%d",i);
1222 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1223 transformation->SetParams(0,0.3,0,&fpar);
1224 kalmanFit->AddCalibration(transformation);
1225 sprintf(tname,"tTPClocalRzIROCOROCSideC_Cos%d",i);
1226 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1227 transformation->SetParams(0,0.3,0,&fpar);
1228 kalmanFit->AddCalibration(transformation);
1231 for (Int_t i=1; i<=nsin;i++){
1233 fpar[0]=0; // cos frequency
1234 fpar[1]=i; // sinus frequency
1235 fpar[2]=1; // relative misalignment
1239 sprintf(tname,"tTPClocalLxIROCOROCSideA_Sin%d",i);
1240 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1241 transformation->SetParams(0,0.03,0,&fpar);
1242 kalmanFit->AddCalibration(transformation);
1243 sprintf(tname,"tTPClocalLxIROCOROCSideC_Sin%d",i);
1244 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1245 transformation->SetParams(0,0.03,0,&fpar);
1246 kalmanFit->AddCalibration(transformation);
1251 sprintf(tname,"tTPClocalLyIROCOROCSideA_Sin%d",i);
1252 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1253 transformation->SetParams(0,0.03,0,&fpar);
1254 kalmanFit->AddCalibration(transformation);
1255 sprintf(tname,"tTPClocalLyIROCOROCSideC_Sin%d",i);
1256 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1257 transformation->SetParams(0,0.03,0,&fpar);
1258 kalmanFit->AddCalibration(transformation);
1264 sprintf(tname,"tTPClocalRzIROCOROCSideA_Sin%d",i);
1265 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1266 transformation->SetParams(0,0.3,0,&fpar);
1267 kalmanFit->AddCalibration(transformation);
1268 sprintf(tname,"tTPClocalRzIROCOROCSideC_Sin%d",i);
1269 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1270 transformation->SetParams(0,0.3,0,&fpar);
1271 kalmanFit->AddCalibration(transformation);
1275 void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1279 AliTPCTransformation * transformation=0;
1281 for (Int_t i=0; i<=ncos;i++){
1283 fpar[0]=i; // cos frequency
1284 fpar[1]=0; // no sinus
1285 fpar[2]=0; // absolute misalignment
1288 // A side is reference
1290 sprintf(tname,"tTPClocalLxSectorSideA_Cos%d",i);
1291 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1292 transformation->SetParams(0,0.03,0,&fpar);
1293 kalmanFit->AddCalibration(transformation);
1298 sprintf(tname,"tTPClocalLySectorSideA_Cos%d",i);
1299 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1300 transformation->SetParams(0,0.03,0,&fpar);
1301 kalmanFit->AddCalibration(transformation);
1309 sprintf(tname,"tTPClocalLxSectorSideC_Cos%d",i);
1310 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1311 transformation->SetParams(0,0.03,0,&fpar);
1312 kalmanFit->AddCalibration(transformation);
1316 sprintf(tname,"tTPClocalLySectorSideC_Cos%d",i);
1317 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1318 transformation->SetParams(0,0.03,0,&fpar);
1319 kalmanFit->AddCalibration(transformation);
1322 // Z rotation - independent
1324 sprintf(tname,"tTPClocalRzSectorSideA_Cos%d",i);
1325 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1326 transformation->SetParams(0,0.3,0,&fpar);
1327 kalmanFit->AddCalibration(transformation);
1328 sprintf(tname,"tTPClocalRzSectorSideC_Cos%d",i);
1329 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1330 transformation->SetParams(0,0.3,0,&fpar);
1331 kalmanFit->AddCalibration(transformation);
1340 for (Int_t i=1; i<=nsin;i++){
1342 fpar[0]=0; // non cos frequency
1343 fpar[1]=i; // sinus frequency
1344 fpar[2]=0; // absolute misalignment
1348 sprintf(tname,"tTPClocalLxSectorSideA_Sin%d",i);
1349 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1350 transformation->SetParams(0,0.03,0,&fpar);
1351 kalmanFit->AddCalibration(transformation);
1352 sprintf(tname,"tTPClocalLxSectorSideC_Sin%d",i);
1353 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1354 transformation->SetParams(0,0.03,0,&fpar);
1355 kalmanFit->AddCalibration(transformation);
1360 sprintf(tname,"tTPClocalLySectorSideA_Sin%d",i);
1361 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1362 transformation->SetParams(0,0.03,0,&fpar);
1363 kalmanFit->AddCalibration(transformation);
1364 sprintf(tname,"tTPClocalLySectorSideC_Sin%d",i);
1365 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1366 transformation->SetParams(0,0.03,0,&fpar);
1367 kalmanFit->AddCalibration(transformation);
1373 sprintf(tname,"tTPClocalRzSectorSideA_Sin%d",i);
1374 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1375 transformation->SetParams(0,0.3,0,&fpar);
1376 kalmanFit->AddCalibration(transformation);
1377 sprintf(tname,"tTPClocalRzSectorSideC_Sin%d",i);
1378 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1379 transformation->SetParams(0,0.3,0,&fpar);
1380 kalmanFit->AddCalibration(transformation);
1385 for (Int_t i=0;i<8;i++){
1386 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1387 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1388 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1389 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1393 void SelectNonPixelA(){
1394 for (Int_t i=0;i<8;i++){
1395 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1396 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1397 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1398 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1402 void SelectNonPixelC(){
1403 for (Int_t i=0;i<8;i++){
1404 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1405 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1406 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1407 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1412 void DumpQA1D( TObjArray &arrayOut){
1415 TF1 fg("fg","gaus");
1416 TMatrixD sideARMS(8,2);
1417 TMatrixD sideCRMS(8,2);
1418 TMatrixD sideACRMS(8,2);
1424 for (Int_t i=0; i<8;i++){
1425 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1426 his->Fit(&fg,"","", -0.15,0.15);
1427 sideARMS(i,0) = fg.GetParameter(2);
1428 his->SetDirectory(0);
1429 his->SetName(Form("Original SideA_%s",his->GetName()));
1430 his->SetTitle(Form("Original SideA_%s",his->GetTitle()));
1431 arrayOut.AddLast(his);
1432 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1433 his->Fit(&fg,"","", -0.15,0.15);
1434 sideARMS(i,1) = fg.GetParameter(2);
1435 his->SetDirectory(0);
1436 his->SetName(Form("Aligned SideA_%s",his->GetName()));
1437 his->SetTitle(Form("Aligned SideA_%s",his->GetTitle()));
1438 arrayOut.AddLast(his);
1444 for (Int_t i=0; i<8;i++){
1445 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1446 his->Fit(&fg,"","", -0.15,0.15);
1447 sideCRMS(i,0) = fg.GetParameter(2);
1448 his->SetDirectory(0);
1449 his->SetName(Form("Original SideC_%s",his->GetName()));
1450 his->SetTitle(Form("Original SideC_%s",his->GetTitle()));
1451 arrayOut.AddLast(his);
1452 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1453 his->Fit(&fg,"","", -0.15,0.15);
1454 sideCRMS(i,1) = fg.GetParameter(2);
1455 his->SetDirectory(0);
1456 his->SetName(Form("Aligned SideC_%s",his->GetName()));
1457 his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1458 arrayOut.AddLast(his);
1464 for (Int_t i=0; i<8;i++){
1465 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1466 his->Fit(&fg,"","", -0.15,0.15);
1467 sideACRMS(i,0) = fg.GetParameter(2);
1468 his->SetDirectory(0);
1469 his->SetName(Form("Original SideAC_%s",his->GetName()));
1470 his->SetTitle(Form("Original SideAC_%s",his->GetTitle()));
1471 arrayOut.AddLast(his);
1472 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1473 his->Fit(&fg,"","", -0.15,0.15);
1474 sideACRMS(i,1) = fg.GetParameter(2);
1475 his->SetDirectory(0);
1476 his->SetName(Form("Aligned SideC_%s",his->GetName()));
1477 his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1478 arrayOut.AddLast(his);
1486 void MakeFolders(TObjArray * arrayOut){
1489 TFolder *folderBase = new TFolder("TPC align","TPC align");
1492 TString maskDelta[8];
1493 TString maskAlign[2]={"Orig","Alig");
1494 for (Int_t i=0; i<8;i++){
1495 maskDelta[i] = kalmanFitNew->fLinearTrackDelta[i]->GetName();
1498 Int_t entries = arrayOut->GetEntries();
1505 void MergeKalman(const char * list = "kalmanFit.list"){
1510 TString currentFile;
1518 printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1519 if (currentFile.Length()==0) continue;
1520 TFile * ffit = TFile::Open(currentFile.Data());
1521 TEntryList * tlist = (TEntryList*) ffit->Get("eventList");
1522 TMatrixD * cMatrix = (TMatrixD*) ffit->Get("cutMarix");
1523 if (tlist&&cMatrix){
1524 printf("Track entries=%d\n",tlist->GetN());
1525 if (cMatrix->Sum()<=0.00000001) {
1526 printf("Problem with track selection\n");
1530 printf("Problem with track selection\n");
1535 AliTPCkalmanFit * fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFit");
1537 if (!kalmanFitNew) {kalmanFitNew= fit; continue;};
1538 kalmanFitNew->Add(fit);
1539 printf("Selected entries=%f\n",fit->fLinearTrackDelta[0]->GetEntries());
1547 currentFile.ReplaceAll("TPC","TPCOrig");
1548 printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1549 if (currentFile.Length()==0) continue;
1550 ffit = TFile::Open(currentFile.Data());
1551 fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFitOrig");
1553 if (!kalmanFitOrig) {kalmanFitOrig= fit; continue;};
1554 kalmanFitOrig->Add(fit);
1560 // save merged results
1562 TFile f("mergeKalmanFit.root","recreate");
1563 kalmanFitNew->Write("kalmanFitNew");
1564 kalmanFitOrig->Write("kalmanFitOrig");
1571 Example - submit alignment as batch jobs
1576 nclustersSkipOffset=0
1583 while [ $ntracksSkipOffset -lt $ntracksSkip ] ; do
1584 nclustersSkipOffset=0;
1585 while [ $nclustersSkipOffset -lt $nclustersSkip ] ; do
1586 echo Tr $ntracksSkipOffset Cl $nclustersSkipOffset;
1587 ver=kalmanDirTrack$ntracksSkipOffset$nclustersSkipOffset
1588 echo New Directory $ver
1591 cp $bDir/align.txt . ;
1592 ln -sf $bDir/kalmanFitApply.root . ;
1593 echo aliroot -q -b "$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C($ntracks,10000,0,$ndump,$ntracksSkip,$ntracksSkipOffset, $nclustersSkip,$nclustersSkipOffset,$isTest)" ;
1594 bsub -q alice-t3_8h -o `pwd`/output.log command aliroot -q -b "$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C($ntracks,10000,0,$ndump,$ntracksSkip,$ntracksSkipOffset, $nclustersSkip,$nclustersSkipOffset,$isTest)" ;
1595 nclustersSkipOffset=$(( $nclustersSkipOffset + 1 ))
1601 ntracksSkipOffset=$(( $ntracksSkipOffset + 1 ))
1602 echo Tr $ntracksSkipOffset;