3 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
4 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
5 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
7 gROOT->LoadMacro("$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C+");
9 AliTPCTransformation::BuildBasicFormulas();
11 AliXRDPROOFtoolkit tool;
12 chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0,50000);
13 chainPoints->Lookup();
15 chainMS = tool.MakeChainRandom("kalmanFit.list","kf",0,50000);
18 chainFP = tool.MakeChainRandom("kalmanFit.list","filter",0,50000);
21 // CalibAlignKalmanFit(40000,1);
22 // kalmanFit0->DumpCorelation(0.8);
23 // TFile f("kalmanfitTPC.root");
37 #include "TLinearFitter.h"
41 #include "TEntryList.h"
43 #include "THnSparse.h"
46 #include "AliSysInfo.h"
47 #include "AliExternalTrackParam.h"
48 #include "TTreeStream.h"
49 #include "AliTrackPointArray.h"
51 #include "AliTPCTransformation.h"
52 #include "AliTPCkalmanFit.h"
53 #include "AliMathBase.h"
54 #include "AliXRDPROOFtoolkit.h"
62 const Float_t krmsYcut = 0.2; // cluster RMS cut in Y - residuals form local fit
63 const Float_t krmsZcut = 0.2; // cluster RMS cut in Z - residuals from local fit
64 const Float_t krmsYcutGlobal= 0.2; // cluster RMS cut in Y - residuals form global fit
65 const Float_t krmsZcutGlobal= 2.0; // cluster RMS cut in Z - residuals from global fit
66 const Float_t kSigmaCut = 5.; // clusters to be removed
67 const Int_t knclCut = 80; // minimal number of clusters
68 const Double_t kArmCut = 50.; // minimal level arm
69 const Double_t kNsigma = 5.; // cut on track level
71 // mult. scattering cuts
73 TCut cutClY("rms09.fElements[0]<0.15");
74 TCut cuttY("rms09.fElements[1]/rms09.fElements[0]<0.9");
75 TCut cutkY("rms09.fElements[2]<0.015");
76 TCut cutClZ("rms09.fElements[3]<0.2");
77 TCut cuttZ("rms09.fElements[4]/rms09.fElements[3]<0.9");
78 TCut cutkZ("rms09.fElements[5]<0.015");
79 TCut cutMS=cutClY+cuttY+cutkY+cutClZ+cuttZ+cutkZ
81 TMatrixD cutMatrix(4*7,2);
82 const Double_t rmsCut09[6]={0.15,0.9,0.015, 0.2, 0.9, 0.015};
89 Int_t toSkipOffset = 0;
90 Int_t toSkipTrack = 2;
91 Int_t toSkipTrackOffset = 0;
92 Int_t isFilterTest = 0;
96 TCut * cSide[4]={0,0,0,0};
97 TCut *cP3[4]={0,0,0,0};
98 TCut *cSP3[4]={0,0,0,0};
99 TCut *cP4[4]={0,0,0,0};
100 TCut *cM4[4]={0,0,0,0};
101 TCut *cA[4]={0,0,0,0};
109 TChain *chainPoints=0;
111 AliTPCkalmanFit * kalmanFitNew=0;
112 AliTPCkalmanFit * kalmanFitOrig=0;
113 AliTPCkalmanFit * kalmanFitApply=0;
115 AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump);
117 AliTPCkalmanFit * SetupFit();
119 void AddFitFieldCage(AliTPCkalmanFit *kalmanFit);
120 void AddPhiScaling(AliTPCkalmanFit *kalmanFit);
121 void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin );
122 void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
123 void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit);
124 void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit);
125 void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
126 void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
127 void AddDrift(AliTPCkalmanFit *kalmanFit);
130 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream);
132 TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream);
133 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset);
135 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump);
137 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){
176 AliTPCTransformation::BuildBasicFormulas();
178 kalmanFitNew = SetupFit();
179 kalmanFitOrig = SetupFit();
180 kalmanFitNew->Init();
181 kalmanFitOrig->Init();
182 return FitPointsLinear(maxTracks,trackDump);
187 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);
209 cSide[0] = new TCut("cutAA","p0In.fP[1]>0&&p1In.fP[1]>0");
210 cSide[1] = new TCut("cutCC","p0In.fP[1]<0&&p1In.fP[1]<0");
211 cSide[2] = new TCut("cutAC","p0In.fP[1]>0&&p1In.fP[1]<0");
212 cSide[3] = new TCut("cutCA","p0In.fP[1]<0&&p1In.fP[1]>0");
214 TH1F * phisP3 = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
215 TH1F * phisSP3 = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
216 TH1F * phisP4 = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
217 TH1F * phisM4 = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
221 TF1 *fg = new TF1("fg","gaus");
222 for (Int_t iter=0; iter<2; iter++){
223 for (Int_t ic=0;ic<4;ic++){
231 chainPoints->Draw("p0In.fP[3]+p1In.fP[3]>>hhisP3",*cA[ic],"goff");
232 phisP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisP3->GetMaximumBin())-0.003,phisP3->GetBinCenter(phisP3->GetMaximumBin())+0.003);
233 cP3[ic]=new TCut(Form("abs(p0In.fP[3]+p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
234 cutMatrix(7*ic+0,0) = fg->GetParameter(1);
235 cutMatrix(7*ic+0,1) = fg->GetParameter(2);
239 chainPoints->Draw("p0Out.fP[3]-p0In.fP[3]>>hhisSP3",*cA[ic],"goff");
240 phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
241 cSP3[ic]=new TCut(Form("abs(p0Out.fP[3]-p0In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
242 cutMatrix(7*ic+1,0) = fg->GetParameter(1);
243 cutMatrix(7*ic+1,1) = fg->GetParameter(2);
245 chainPoints->Draw("p1Out.fP[3]-p1In.fP[3]>>hhisSP3",*cA[ic],"goff");
246 phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
247 *cSP3[ic]+=Form("abs(p1Out.fP[3]-p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
248 cutMatrix(7*ic+2,0) = fg->GetParameter(1);
249 cutMatrix(7*ic+2,1) = fg->GetParameter(2);
253 chainPoints->Draw("p0Out.fP[4]>>hhisP4",*cA[ic],"goff");
254 phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
255 cP4[ic]=new TCut(Form("abs(p0Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
256 cutMatrix(7*ic+3,0) = fg->GetParameter(1);
257 cutMatrix(7*ic+3,1) = fg->GetParameter(2);
258 chainPoints->Draw("p1Out.fP[4]>>hhisP4",*cA[ic],"goff");
259 phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
260 *cP4[ic]+=Form("abs(p1Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
261 cutMatrix(7*ic+4,0) = fg->GetParameter(1);
262 cutMatrix(7*ic+4,1) = fg->GetParameter(2);
267 chainPoints->Draw("p0Out.fP[4]-p0In.fP[4]>>hhisM4",*cA[ic],"goff");
268 phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
269 cM4[ic]=new TCut(Form("abs(p0Out.fP[4]-p0In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
270 cutMatrix(7*ic+5,0) = fg->GetParameter(1);
271 cutMatrix(7*ic+5,1) = fg->GetParameter(2);
273 chainPoints->Draw("p1Out.fP[4]-p1In.fP[4]>>hhisM4",*cA[ic],"goff");
274 phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
275 *cM4[ic]+=Form("abs(p1Out.fP[4]-p1In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
276 cutMatrix(7*ic+6,0) = fg->GetParameter(1);
277 cutMatrix(7*ic+6,1) = fg->GetParameter(2);
282 *cA[ic]= *cSide[ic]+*cP3[ic]+*cSP3[ic]+*cP4[ic]+*cM4[ic];
287 cutAll = (*cA[0])||(*cA[1])||(*cA[2])||(*cA[3])+"abs(mag)<0.01&&ncont>0&&p.fNPoints>120";
289 delete phisP3; // = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
290 delete phisSP3; // = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
291 delete phisP4;// = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
292 delete phisM4;// = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
296 chainPoints->Draw(">>listEL",cutAll,"entryList");
297 elist = (TEntryList*)gDirectory->Get("listEL");
298 chainPoints->SetEntryList(elist);
299 elist->SetDirectory(0);
304 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump){
308 // create debug streeamers
309 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPC.root");
310 TTreeSRedirector *pcstreamOrig = new TTreeSRedirector("kalmanfitTPCOrig.root");
311 pcstream->GetFile()->cd();
312 cutMatrix.Write("cutMarix");
313 elist->Write("eventList");
316 AliTrackPointArray *pointsNS=0;
319 AliExternalTrackParam *param0=0;
320 AliExternalTrackParam *param1=0;
321 chainPoints->SetBranchAddress("p.",&pointsNS);
322 chainPoints->SetBranchAddress("p0In.",¶m0);
323 chainPoints->SetBranchAddress("p1In.",¶m1);
324 chainPoints->SetBranchAddress("mag",&mag);
325 chainPoints->SetBranchAddress("time",&time);
327 printf("\n*\n*\n*Selected entries = %d\n*\n*\n*",Int_t(elist->GetN()));
330 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
331 if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
332 Int_t entry=chainPoints->GetEntryNumber(itrack);
333 chainPoints->GetEntry(entry);
334 if (accepted>maxTracks) break;
336 AliTrackPointArray *points = AliTPCkalmanFit::SortPoints(*pointsNS);
337 if (kalmanFitApply) kalmanFitApply->ApplyCalibration(points,-1.);
339 // estimate and filter scattering
341 TVectorD *vecRMS09 = EstimateScatteringKalmanLinear(*points,*param0,*param1,pcstream);
342 if (!vecRMS09) continue;
344 if ((*vecRMS09)[0] >rmsCut09[0]) isOK=kFALSE;
345 if ((*vecRMS09)[1]/(*vecRMS09)[0] >rmsCut09[1]) isOK=kFALSE;;
346 if ((*vecRMS09)[2] >rmsCut09[2]) isOK=kFALSE;
347 if ((*vecRMS09)[3] >rmsCut09[3]) isOK=kFALSE;
348 if ((*vecRMS09)[4]/(*vecRMS09)[0] >rmsCut09[4]) isOK=kFALSE;
349 if ((*vecRMS09)[5] >rmsCut09[5]) isOK=kFALSE;
350 if (!isOK || isFilterTest) {
354 kalmanFitNew->PropagateTime(time);
357 for (Int_t idir=-1; idir<=1; idir++){
358 AliTrackPointArray *pointsF = FilterPoints(*points,idir, pcstream);
359 if (!pointsF) continue;
360 AliTrackPointArray *spointsF = 0;
361 // we skip points for alignemnt but not for QA
362 if (idir==0) spointsF = SkipPoints(*pointsF, toSkip*2, toSkipOffset);
363 if (idir!=0) spointsF = SkipPoints(*pointsF, toSkip, toSkipOffset);
364 if (idir==0) accepted++;
366 if (accepted%50==0) {
367 kalmanFitNew->FitTrackLinear(*pointsF, pcstream);
369 if (idir==0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
370 if (idir!=0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
372 if (idir==0) kalmanFitNew->DumpTrackLinear(*pointsF,pcstream);
373 if (idir==0) kalmanFitOrig->DumpTrackLinear(*pointsF,pcstreamOrig);
374 if (accepted%trackDump==0) {
375 printf("%d\n", accepted);
377 AliSysInfo::AddStamp("trackFit", accepted,itrack);
383 pcstream->GetFile()->cd();
384 kalmanFitNew->Write("kalmanFit");
385 pcstreamOrig->GetFile()->cd();
386 kalmanFitOrig->Write("kalmanFitOrig");
387 pcstreamOrig->GetFile()->cd();
388 if (kalmanFitApply) kalmanFitApply->Write("kalmanFitApply");
395 void QAPointsLinear(Int_t maxTracks, Int_t trackDump){
397 // check the consistency of kalman fit
398 // Apply transformation
400 // create debug streeamers
401 TTreeSRedirector *pcstreamNonCalib = new TTreeSRedirector("kalmanfitTPCQANonCalib.root");
402 TTreeSRedirector *pcstreamCalib = new TTreeSRedirector("kalmanfitTPCQACalib.root");
403 TTreeSRedirector *pcstream=0;
404 AliTPCkalmanFit *kalmanFitters[6]={0,0,0,0,0,0};
405 for (Int_t i=0;i<6;i++){
406 kalmanFitters[i]=SetupFit();
410 AliTrackPointArray *points=0;
411 AliExternalTrackParam *param0=0;
412 AliExternalTrackParam *param1=0;
415 chainPoints->SetBranchAddress("p.",&points);
416 chainPoints->SetBranchAddress("mag",&mag);
417 chainPoints->SetBranchAddress("time",&time);
418 chainPoints->SetBranchAddress("p0In.",¶m0);
419 chainPoints->SetBranchAddress("p1In.",¶m1);
425 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
426 if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
427 Int_t entry=chainPoints->GetEntryNumber(itrack);
428 chainPoints->GetEntry(entry);
429 if (accepted>maxTracks) break;
432 AliTrackPointArray pointsCalib(*points);
433 for (Int_t iscalib=0; iscalib<1;iscalib++){
434 if (iscalib>0) kalmanFitNew->ApplyCalibration(&pointsCalib,-1.);
435 if (iscalib==0) pcstream=pcstreamNonCalib;
436 if (iscalib>0) pcstream=pcstreamCalib;
437 for (Int_t idir=-1; idir<=1; idir++){
438 AliTrackPointArray *pointsF = FilterPoints(pointsCalib,idir, pcstream);
439 if (!pointsF) continue;
441 if (idir==0) accepted++;
442 kalmanFitters[iscalib*3+idir+1]->DumpTrackLinear(*pointsF,0);
443 EstimateScatteringKalmanLinear(*pointsF,*param0,*param1,pcstream);
447 if (accepted%trackDump==0) {
448 printf("%d\n", accepted);
451 pcstreamCalib->GetFile()->cd();
452 kalmanFitters[0]->Write("fitUpNonCalib");
453 kalmanFitters[1]->Write("fitUpDownNonCalib");
454 kalmanFitters[2]->Write("fitDownNonCalib");
455 kalmanFitters[3]->Write("fitUpCalib");
456 kalmanFitters[4]->Write("fitUpDownCalib");
457 kalmanFitters[5]->Write("fitDownCalib");
458 delete pcstreamCalib;
459 delete pcstreamNonCalib;
463 void TestScattering(Int_t maxTracks, Int_t trackDump){
465 // test Multiple scattering algorithm
466 // Apply transformation
468 // create debug streeamers
469 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPCMS.root");
472 AliTrackPointArray *points=0;
473 AliExternalTrackParam *param0=0;
474 AliExternalTrackParam *param1=0;
477 chainPoints->SetBranchAddress("p.",&points);
478 chainPoints->SetBranchAddress("mag",&mag);
479 chainPoints->SetBranchAddress("time",&time);
480 chainPoints->SetBranchAddress("p0In.",¶m0);
481 chainPoints->SetBranchAddress("p1In.",¶m1);
484 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
485 Int_t entry=chainPoints->GetEntryNumber(itrack);
486 chainPoints->GetEntry(entry);
487 if (accepted>maxTracks) break;
489 AliTrackPointArray *pointsSorted = AliTPCkalmanFit::SortPoints(*points);
490 EstimateScatteringKalmanLinear(*pointsSorted,*param0,*param1,pcstream);
492 if (accepted%trackDump==0) {
493 printf("%d\n", accepted);
504 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset){
506 // create new array with skipped points
508 Int_t npoints = points.GetNPoints();
509 Int_t npointsF = (npoints-nskipOffset-1)/nskip;
511 AliTrackPointArray *pointsF= new AliTrackPointArray(npointsF);
513 for (Int_t ipoint=nskipOffset; ipoint<npoints; ipoint+=nskip){
515 if (!points.GetPoint(point,ipoint)) continue;
516 pointsF->AddPoint(used,&point);
518 if (used==npointsF) break;
525 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream){
527 // Filter points - input points for KalmanFilter
530 TLinearFitter lfitY(2,"pol1");
531 TLinearFitter lfitZ(2,"pol1");
535 lfitY.StoreData(kTRUE);
536 lfitZ.StoreData(kTRUE);
537 Int_t npoints = points.GetNPoints();
538 if (npoints<2) return 0;
539 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
540 Double_t ca = TMath::Cos(currentAlpha);
541 Double_t sa = TMath::Sin(currentAlpha);
543 // 1.b Fit the track in the rotated frame - MakeSeed
545 Double_t maxX =-10000, minX=10000;
546 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
547 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
548 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
549 Double_t rz = points.GetZ()[ipoint];
550 if (dir== 1 && rx<0) continue;
551 if (dir==-1 && rx>0) continue;
552 if (maxX<rx) maxX=rx;
553 if (minX>rx) minX=rx;
554 lfitY.AddPoint(&rx,ry,1);
555 lfitZ.AddPoint(&rx,rz,1);
557 if (TMath::Abs(maxX-minX)<kArmCut) return 0;
558 if (lfitY.GetNpoints()<knclCut) return 0;
562 lfitY.GetParameters(vecY);
563 lfitZ.GetParameters(vecZ);
565 Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
566 Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
567 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
568 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
574 AliTrackPointArray *pointsF=0;
575 for (Int_t iter=0; iter<2;iter++){
576 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
578 if (!points.GetPoint(point,ipoint)) continue;
579 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
580 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
581 Double_t rz = points.GetZ()[ipoint];
582 if (dir== 1 && rx<0) continue;
583 if (dir==-1 && rx>0) continue;
584 Double_t erry = TMath::Sqrt(chi2Y);
585 Double_t errz = TMath::Sqrt(chi2Z);
586 Double_t fy = vecY[0]+vecY[1]*rx;
587 Double_t fz = vecZ[0]+vecZ[1]*rx;
588 if (TMath::Abs(fy-ry)>erry*kSigmaCut) continue;
589 if (TMath::Abs(fz-rz)>errz*kSigmaCut) continue;
591 if (pointsF) pointsF->AddPoint(toBeUsed,&point);
595 (*pcstream)<<"filter"<<
597 "accepted="<<accepted<<
606 if (accepted<knclCut) break;
607 if (iter==0) pointsF = new AliTrackPointArray(toBeUsed);
615 AliTrackPointArray * SortPoints(AliTrackPointArray &points){
617 //Creates the array - points sorted according radius - neccessay for kalman fit
620 // 0. choose the frame - rotation angle
622 Int_t npoints = points.GetNPoints();
623 if (npoints<1) return 0;
624 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
625 Double_t ca = TMath::Cos(currentAlpha);
626 Double_t sa = TMath::Sin(currentAlpha);
628 // 1. sort the points
630 Double_t *rxvector = new Double_t[npoints];
631 Int_t *indexes = new Int_t[npoints];
632 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
633 rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
635 TMath::Sort(npoints, rxvector,indexes,kFALSE);
637 AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
638 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
639 if (!points.GetPoint(point,indexes[ipoint])) continue;
640 pointsSorted->AddPoint(ipoint,&point);
647 TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream){
649 // Algorithm - 0. Fit the track forward and backward
650 // - 1. Store the current parameters in each point
652 const Int_t kMinPoints= 70;
653 const Double_t kResY = 0.1;
654 const Double_t kResZ = 0.1;
655 const Double_t kMisY = 0.02;
656 const Double_t kMisZ = 0.02;
657 const Double_t kLArm = 120.;
658 const Double_t kACsideFac = 10.;
659 const Double_t kMS0 = kResY/20.;
661 Int_t npoints = points.GetNPoints();
662 if (npoints<kMinPoints) return 0;
664 TVectorD *vecPos[11]={0,0,0,0,0,0,0,0,0,0,0};
665 for (Int_t i=0;i<11;i++){
666 vecPos[i]=new TVectorD(npoints);
670 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
671 Double_t ca = TMath::Cos(currentAlpha);
672 Double_t sa = TMath::Sin(currentAlpha);
675 TMatrixD trParamY(2,1),trCovarY(2,2);
676 TMatrixD trParamZ(2,1),trCovarZ(2,2);
680 TMatrixD matHk(kNmeas,nelem); // vector to mesurement
681 TMatrixD matHkT(nelem,kNmeas); // helper matrix Hk transpose
682 TMatrixD matFk(nelem,nelem);
683 TMatrixD matFkT(nelem,nelem);
684 TMatrixD vecYk(kNmeas,1); // Innovation or measurement residual
685 TMatrixD vecZk(kNmeas,1); // Innovation or measurement residual
686 TMatrixD measR(kNmeas,kNmeas);
687 TMatrixD matSk(kNmeas,kNmeas); // Innovation (or residual) covariance
688 TMatrixD matKk(nelem,kNmeas); // Optimal Kalman gain
689 TMatrixD covXk2(nelem,nelem); // helper matrix
690 TMatrixD covXk3(nelem,nelem); // helper matrix
691 TMatrixD mat1(nelem,nelem);
692 mat1(0,0)=1.; mat1(0,1)=0.;
693 mat1(1,0)=0.; mat1(1,1)=1.;
699 for (Int_t idir=0; idir<2;idir++){
702 for (Int_t ip=0; ip<npoints; ip++){
704 if (idir>0) ipoint = npoints-ip-1;
705 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
706 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
707 Double_t rz = points.GetZ()[ipoint];
708 Int_t volId= points.GetVolumeID()[ipoint];
711 // set initital parameters and covariance - use first and middle point
712 Double_t rxm = ca*points.GetX()[npoints/2]+sa*points.GetY()[npoints/2];
713 Double_t rym = -sa*points.GetX()[npoints/2]+ca*points.GetY()[npoints/2];
714 Double_t rzm = points.GetZ()[npoints/2];
716 trParamY(1,0) = (rym-ry)/(rxm-rx);
718 trParamZ(1,0) = (rzm-rz)/(rxm-rx);
720 trCovarY(0,0) = kResY*kResY;
721 trCovarY(1,1) = (kResY*kResY)/((rxm-rx)*(rxm-rx));
722 trCovarZ(0,0) = kResZ*kResZ;
723 trCovarZ(1,1) = (kResZ*kResZ)/((rxm-rx)*(rxm-rx));
730 if ((volId%36)<18 && (lastVolId%36)>=18){
732 trCovarY(0,0)+=kMisY*kMisY*kACsideFac;
733 trCovarZ(0,0)+=kMisZ*kMisZ*kACsideFac;
734 trCovarY(1,1)+=kACsideFac*(kMisY*kMisY)/(kLArm*kLArm);;
735 trCovarZ(1,1)+=kACsideFac*(kMisZ*kMisZ)/(kLArm*kLArm);
738 if (volId!=lastVolId){
740 trCovarY(0,0)+=kMisY*kMisY;
741 trCovarZ(0,0)+=kMisZ*kMisZ;
742 trCovarY(1,1)+=(kMisY*kMisY)/(kLArm*kLArm);;
743 trCovarZ(1,1)+=(kMisZ*kMisZ)/(kLArm*kLArm);
749 Double_t deltaX=rx-lastX;
750 trParamY(0,0)+=deltaX*trParamY(1,0);
751 trParamZ(0,0)+=deltaX*trParamZ(1,0);
752 matFk(0,0)=1; matFk(0,1)=deltaX;
753 matFk(1,0)=0; matFk(1,1)=1;
754 matFkT=matFk.T(); matFk.T();
755 covXk2=matFk*trCovarY*matFkT;
757 covXk2=matFk*trCovarZ*matFkT;
760 // multiple scattering
761 trCovarY(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
762 trCovarZ(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
767 for (Int_t coord=0; coord<2;coord++){
768 TMatrixD* pvecXk = (coord==0)? &trParamY: &trParamZ;
769 TMatrixD* pcovXk = (coord==0)? &trCovarY: &trCovarZ;
771 TMatrixD& vecXk = *pvecXk;
772 TMatrixD& covXk = *pcovXk;
773 measR(0,0) = (coord==0) ? kResY:kResZ;
774 vecZk(0,0) = (coord==0) ? ry:rz;
776 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
777 matHkT=matHk.T(); matHk.T();
778 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
780 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
782 covXk2= (mat1-(matKk*matHk));
783 covXk3 = covXk2*covXk;
784 vecXk += matKk*vecYk; // updated vector
790 (*vecPos[0])[ipoint]=rx;
791 (*vecPos[1])[ipoint]=ry;
792 (*vecPos[2])[ipoint]=rz;
794 (*vecPos[4*idir+0+3])[ipoint]=trParamY(0,0);
795 (*vecPos[4*idir+1+3])[ipoint]=trParamY(1,0);
797 (*vecPos[4*idir+2+3])[ipoint]=trParamZ(0,0);
798 (*vecPos[4*idir+3+3])[ipoint]=trParamZ(1,0);
804 TVectorD vec(npoints);
806 TVectorD rms09(6); // robust RMS - fraction 0.9
807 TVectorD mean09(6); // robust RMS - fraction 0.9
809 Int_t npoints09 = Int_t(npoints*0.9);
811 vec=(*(vecPos[3])-*(vecPos[1]));
812 rms[0]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms cluster y
813 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
816 vec=(*(vecPos[7])-*(vecPos[3]));
817 rms[1]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track y
818 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
821 vec=(*(vecPos[8])-*(vecPos[4]));
822 rms[2]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track ky
823 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
827 vec=(*(vecPos[5])-*(vecPos[2]));
828 rms[3]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms cluster z
829 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
832 vec=(*(vecPos[9])-*(vecPos[5]));
833 rms[4]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track z
834 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
837 vec=(*(vecPos[10])-*(vecPos[6]));
838 rms[5]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track kz
839 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
850 "mean09.="<<&mean09<<
856 for (Int_t i=0;i<11;i++){
859 return new TVectorD(rms09);
870 void AddFitFieldCage(AliTPCkalmanFit *kalmanFit){
872 // Add radial scaling due field cage
875 AliTPCTransformation * transformation=0;
878 // linear R scaling and shift
880 for (Int_t iside=0; iside<=1; iside++)
881 for (Int_t ipolR=0; ipolR<2; ipolR++){
882 for (Int_t ipolZ=0; ipolZ<3; ipolZ++){
885 if (ipolR+ipolZ==0) continue;
886 sprintf(tname,"tTPCscalingRPolR%dDr%dSide%d",ipolR,ipolZ,iside);
887 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0, 1);
888 transformation->SetParams(0,0.2,0,&fpar);
889 kalmanFit->AddCalibration(transformation);
896 for (Int_t iside=0; iside<=1; iside++)
897 for (Int_t ipol=0; ipol<3; ipol++){
899 sprintf(tname,"tTPCscalingRIFC%dSide%d",ipol,iside);
900 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRIFC",0,0, 1);
901 transformation->SetParams(0,0.2,0,&fpar);
902 kalmanFit->AddCalibration(transformation);
907 for (Int_t iside=0; iside<=1; iside++)
908 for (Int_t ipol=0; ipol<3; ipol++){
911 sprintf(tname,"tTPCscalingROFC%dSide%d",ipol,iside);
912 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingROFC",0,0, 1);
913 transformation->SetParams(0,0.2,0,&fpar);
914 kalmanFit->AddCalibration(transformation);
919 void AddPhiScaling(AliTPCkalmanFit *kalmanFit){
921 // Add linear local phi scaling -
922 // separate IROC/OROC - A side/C side
923 // "tscalingLocalPhiIROC"
924 // "tscalingLocalPhiOROC"
928 for (Int_t i=0;i<72;i++){
937 AliTPCTransformation * transformation=0;
939 transformation = new AliTPCTransformation("tscalingLocalPhiIROC", new TBits(maskInner), 0,"TPCscalingPhiLocal",0, 1);
940 transformation->SetParams(0,0.02,0,&fpar);
941 kalmanFit->AddCalibration(transformation);
942 transformation = new AliTPCTransformation("tscalingLocalPhiOROC", new TBits(maskOuter), 0,"TPCscalingPhiLocal",0, 1);
943 transformation->SetParams(0,0.02,0,&fpar);
944 kalmanFit->AddCalibration(transformation);
948 void AddDrift(AliTPCkalmanFit *kalmanFit){
950 // Add drift velocity transformation
953 AliTPCTransformation * transformation=0;
955 transformation = new AliTPCTransformation("tTPCscalingZDrift", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDr", 0);
956 transformation->SetParams(0,4.,1.,&fpar);
957 kalmanFit->AddCalibration(transformation);
959 transformation = new AliTPCTransformation("tTPCscalingZDriftGy", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDrGy", 0);
960 transformation->SetParams(0,0.2,0.0,&fpar);
961 kalmanFit->AddCalibration(transformation);
967 void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
972 fpar[0]=0; fpar[1]=0; fpar[2]=0;
974 AliTPCTransformation * transformation=0;
978 for (Int_t i=0; i<=ncos;i++){
979 fpar[0]=i; // cos frequency
980 fpar[1]=0; // no sinus
981 fpar[2]=1; // relative misalignment
983 sprintf(tname,"tTPCDeltaZIROCOROCSideA_Cos%d",i);
984 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
985 transformation->SetParams(0,0.03,0,&fpar);
986 kalmanFit->AddCalibration(transformation);
988 sprintf(tname,"tTPCDeltaZIROCOROCSideC_Cos%d",i);
989 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
990 transformation->SetParams(0,0.03,0,&fpar);
991 kalmanFit->AddCalibration(transformation);
995 fpar[2]=0; // absolute misalignment
996 sprintf(tname,"tTPCDeltaZSectorSideA_Cos%d",i);
997 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
998 transformation->SetParams(0,0.1,0,&fpar);
999 kalmanFit->AddCalibration(transformation);
1001 sprintf(tname,"tTPCDeltaZSectorSideC_Cos%d",i);
1002 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1003 transformation->SetParams(0,0.1,0,&fpar);
1004 kalmanFit->AddCalibration(transformation);
1007 for (Int_t i=1; i<=nsin;i++){
1008 fpar[0]=0; // cos frequency
1009 fpar[1]=i; // sinus frequncy
1010 fpar[2]=1; // relative misalignment
1012 sprintf(tname,"tTPCDeltaZIROCOROCSideA_Sin%d",i);
1013 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
1014 transformation->SetParams(0,0.03,0,&fpar);
1015 kalmanFit->AddCalibration(transformation);
1017 sprintf(tname,"tTPCDeltaZIROCOROCSideC_Sin%d",i);
1018 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1019 transformation->SetParams(0,0.03,0,&fpar);
1020 kalmanFit->AddCalibration(transformation);
1024 fpar[2]=0; // absolute misalignment
1025 sprintf(tname,"tTPCDeltaZSectorSideA_Sin%d",i);
1026 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
1027 transformation->SetParams(0,0.1,0,&fpar);
1028 kalmanFit->AddCalibration(transformation);
1030 sprintf(tname,"tTPCDeltaZSectorSideC_Sin%d",i);
1031 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1032 transformation->SetParams(0,0.1,0,&fpar);
1033 kalmanFit->AddCalibration(transformation);
1041 void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1043 // z tilting absolute (sector) and relative (IROC-OROC)
1047 fpar[0]=0; fpar[1]=0; fpar[2]=0;
1049 AliTPCTransformation * transformation=0;
1053 for (Int_t i=0; i<=ncos;i++){
1054 fpar[0]=i; // cos frequency
1055 fpar[1]=0; // sinus frequency
1056 fpar[2]=1; // relative misalignment
1058 sprintf(tname,"tTPCTiltingZIROCOROCSideA_Cos%d",i);
1059 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1060 transformation->SetParams(0,0.03,0,&fpar);
1061 kalmanFit->AddCalibration(transformation);
1063 sprintf(tname,"tTPCTiltingZIROCOROCSideC_Cos%d",i);
1064 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1065 transformation->SetParams(0,0.03,0,&fpar);
1066 kalmanFit->AddCalibration(transformation);
1070 fpar[2]=0; // absolute misalignment
1071 sprintf(tname,"tTPCTiltingZSectorSideA_Cos%d",i);
1072 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1073 transformation->SetParams(0,0.1,0,&fpar);
1074 kalmanFit->AddCalibration(transformation);
1076 sprintf(tname,"tTPCTiltingZSectorSideC_Cos%d",i);
1077 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1078 transformation->SetParams(0,0.1,0,&fpar);
1079 kalmanFit->AddCalibration(transformation);
1082 for (Int_t i=1; i<=nsin;i++){
1083 fpar[0]=0; // cos frequency
1084 fpar[1]=i; // sinus frequncy
1085 fpar[2]=1; // relative misalignment
1087 sprintf(tname,"tTPCTiltingZIROCOROCSideA_Sin%d",i);
1088 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1089 transformation->SetParams(0,0.03,0,&fpar);
1090 kalmanFit->AddCalibration(transformation);
1092 sprintf(tname,"tTPCTiltingZIROCOROCSideC_Sin%d",i);
1093 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1094 transformation->SetParams(0,0.03,0,&fpar);
1095 kalmanFit->AddCalibration(transformation);
1099 fpar[2]=0; // absolute misalignment
1100 sprintf(tname,"tTPCTiltingZSectorSideA_Sin%d",i);
1101 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1102 transformation->SetParams(0,0.1,0,&fpar);
1103 kalmanFit->AddCalibration(transformation);
1105 sprintf(tname,"tTPCTiltingZSectorSideC_Sin%d",i);
1106 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1107 transformation->SetParams(0,0.1,0,&fpar);
1108 kalmanFit->AddCalibration(transformation);
1114 void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit){
1119 AliTPCTransformation * transformation=0;
1120 TBits maskInnerA(72);
1121 TBits maskInnerC(72);
1122 for (Int_t i=0;i<72;i++){
1124 if (i%36<18) maskInnerA[i]=kTRUE;
1125 if (i%36>=18) maskInnerC[i]=kTRUE;
1130 transformation = new AliTPCTransformation("tTPCDeltaLxIROCA", new TBits(maskInnerA), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1131 transformation->SetParams(0,0.2,0,&fpar);
1132 kalmanFit->AddCalibration(transformation);
1133 transformation = new AliTPCTransformation("tTPCDeltaLxIROCC", new TBits(maskInnerC), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1134 transformation->SetParams(0,0.2,0,&fpar);
1135 kalmanFit->AddCalibration(transformation);
1137 transformation = new AliTPCTransformation("tTPCDeltaLyIROCA", new TBits(maskInnerA), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1138 transformation->SetParams(0,0.2,0,&fpar);
1139 kalmanFit->AddCalibration(transformation);
1140 transformation = new AliTPCTransformation("tTPCDeltaLyIROCC", new TBits(maskInnerC), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1141 transformation->SetParams(0,0.2,0,&fpar);
1142 kalmanFit->AddCalibration(transformation);
1145 void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit){
1150 AliTPCTransformation * transformation=0;
1153 for (Int_t isec=0; isec<36;isec++){
1158 sprintf(tname,"tTPClocalLxIROC%d",isec);
1159 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1160 transformation->SetParams(0,0.2,0,&fpar);
1161 kalmanFit->AddCalibration(transformation);
1163 sprintf(tname,"tTPClocalLyIROC%d",isec);
1164 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1165 transformation->SetParams(0,0.2,0,&fpar);
1166 kalmanFit->AddCalibration(transformation);
1168 sprintf(tname,"tTPClocalRzIROC%d",isec);
1169 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0, 0);
1170 transformation->SetParams(0,0.2,0,&fpar);
1171 kalmanFit->AddCalibration(transformation);
1175 for (Int_t isec=0; isec<36;isec++){
1176 if (isec%18==fixSector) continue;
1179 mask[isec+36]=kTRUE;
1182 sprintf(tname,"tTPClocalLxSector%d",isec);
1183 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1184 transformation->SetParams(0,0.2,0,&fpar);
1185 kalmanFit->AddCalibration(transformation);
1187 sprintf(tname,"tTPClocalLySector%d",isec);
1188 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1189 transformation->SetParams(0,0.2,0,&fpar);
1190 kalmanFit->AddCalibration(transformation);
1192 sprintf(tname,"tTPClocalRzSector%d",isec);
1193 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0, 0);
1194 transformation->SetParams(0,0.2,0,&fpar);
1195 kalmanFit->AddCalibration(transformation);
1202 void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1207 AliTPCTransformation * transformation=0;
1209 for (Int_t i=0; i<=ncos;i++){
1211 fpar[0]=i; // cos frequency
1212 fpar[1]=0; // no sinus
1213 fpar[2]=1; // relative misalignment
1217 sprintf(tname,"tTPClocalLxIROCOROCSideA_Cos%d",i);
1218 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1219 transformation->SetParams(0,0.03,0,&fpar);
1220 kalmanFit->AddCalibration(transformation);
1221 sprintf(tname,"tTPClocalLxIROCOROCSideC_Cos%d",i);
1222 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1223 transformation->SetParams(0,0.03,0,&fpar);
1224 kalmanFit->AddCalibration(transformation);
1229 sprintf(tname,"tTPClocalLyIROCOROCSideA_Cos%d",i);
1230 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1231 transformation->SetParams(0,0.03,0,&fpar);
1232 kalmanFit->AddCalibration(transformation);
1233 sprintf(tname,"tTPClocalLyIROCOROCSideC_Cos%d",i);
1234 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1235 transformation->SetParams(0,0.03,0,&fpar);
1236 kalmanFit->AddCalibration(transformation);
1242 sprintf(tname,"tTPClocalRzIROCOROCSideA_Cos%d",i);
1243 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1244 transformation->SetParams(0,0.3,0,&fpar);
1245 kalmanFit->AddCalibration(transformation);
1246 sprintf(tname,"tTPClocalRzIROCOROCSideC_Cos%d",i);
1247 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1248 transformation->SetParams(0,0.3,0,&fpar);
1249 kalmanFit->AddCalibration(transformation);
1252 for (Int_t i=1; i<=nsin;i++){
1254 fpar[0]=0; // cos frequency
1255 fpar[1]=i; // sinus frequency
1256 fpar[2]=1; // relative misalignment
1260 sprintf(tname,"tTPClocalLxIROCOROCSideA_Sin%d",i);
1261 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1262 transformation->SetParams(0,0.03,0,&fpar);
1263 kalmanFit->AddCalibration(transformation);
1264 sprintf(tname,"tTPClocalLxIROCOROCSideC_Sin%d",i);
1265 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1266 transformation->SetParams(0,0.03,0,&fpar);
1267 kalmanFit->AddCalibration(transformation);
1272 sprintf(tname,"tTPClocalLyIROCOROCSideA_Sin%d",i);
1273 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1274 transformation->SetParams(0,0.03,0,&fpar);
1275 kalmanFit->AddCalibration(transformation);
1276 sprintf(tname,"tTPClocalLyIROCOROCSideC_Sin%d",i);
1277 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1278 transformation->SetParams(0,0.03,0,&fpar);
1279 kalmanFit->AddCalibration(transformation);
1285 sprintf(tname,"tTPClocalRzIROCOROCSideA_Sin%d",i);
1286 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1287 transformation->SetParams(0,0.3,0,&fpar);
1288 kalmanFit->AddCalibration(transformation);
1289 sprintf(tname,"tTPClocalRzIROCOROCSideC_Sin%d",i);
1290 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1291 transformation->SetParams(0,0.3,0,&fpar);
1292 kalmanFit->AddCalibration(transformation);
1296 void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1301 AliTPCTransformation * transformation=0;
1303 for (Int_t i=0; i<=ncos;i++){
1305 fpar[0]=i; // cos frequency
1306 fpar[1]=0; // no sinus
1307 fpar[2]=0; // absolute misalignment
1310 // A side is reference
1312 sprintf(tname,"tTPClocalLxSectorSideA_Cos%d",i);
1313 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1314 transformation->SetParams(0,0.03,0,&fpar);
1315 kalmanFit->AddCalibration(transformation);
1320 sprintf(tname,"tTPClocalLySectorSideA_Cos%d",i);
1321 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1322 transformation->SetParams(0,0.03,0,&fpar);
1323 kalmanFit->AddCalibration(transformation);
1331 sprintf(tname,"tTPClocalLxSectorSideC_Cos%d",i);
1332 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1333 transformation->SetParams(0,0.03,0,&fpar);
1334 kalmanFit->AddCalibration(transformation);
1338 sprintf(tname,"tTPClocalLySectorSideC_Cos%d",i);
1339 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1340 transformation->SetParams(0,0.03,0,&fpar);
1341 kalmanFit->AddCalibration(transformation);
1344 // Z rotation - independent
1346 sprintf(tname,"tTPClocalRzSectorSideA_Cos%d",i);
1347 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1348 transformation->SetParams(0,0.3,0,&fpar);
1349 kalmanFit->AddCalibration(transformation);
1350 sprintf(tname,"tTPClocalRzSectorSideC_Cos%d",i);
1351 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1352 transformation->SetParams(0,0.3,0,&fpar);
1353 kalmanFit->AddCalibration(transformation);
1362 for (Int_t i=1; i<=nsin;i++){
1364 fpar[0]=0; // non cos frequency
1365 fpar[1]=i; // sinus frequency
1366 fpar[2]=0; // absolute misalignment
1370 sprintf(tname,"tTPClocalLxSectorSideA_Sin%d",i);
1371 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1372 transformation->SetParams(0,0.03,0,&fpar);
1373 kalmanFit->AddCalibration(transformation);
1374 sprintf(tname,"tTPClocalLxSectorSideC_Sin%d",i);
1375 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1376 transformation->SetParams(0,0.03,0,&fpar);
1377 kalmanFit->AddCalibration(transformation);
1382 sprintf(tname,"tTPClocalLySectorSideA_Sin%d",i);
1383 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1384 transformation->SetParams(0,0.03,0,&fpar);
1385 kalmanFit->AddCalibration(transformation);
1386 sprintf(tname,"tTPClocalLySectorSideC_Sin%d",i);
1387 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1388 transformation->SetParams(0,0.03,0,&fpar);
1389 kalmanFit->AddCalibration(transformation);
1395 sprintf(tname,"tTPClocalRzSectorSideA_Sin%d",i);
1396 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1397 transformation->SetParams(0,0.3,0,&fpar);
1398 kalmanFit->AddCalibration(transformation);
1399 sprintf(tname,"tTPClocalRzSectorSideC_Sin%d",i);
1400 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1401 transformation->SetParams(0,0.3,0,&fpar);
1402 kalmanFit->AddCalibration(transformation);
1407 for (Int_t i=0;i<8;i++){
1408 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1409 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1410 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1411 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1415 void SelectNonPixelA(){
1416 for (Int_t i=0;i<8;i++){
1417 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1418 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1419 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1420 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1424 void SelectNonPixelC(){
1425 for (Int_t i=0;i<8;i++){
1426 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1427 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1428 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1429 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1434 void DumpQA1D( TObjArray &arrayOut){
1438 TF1 fg("fg","gaus");
1439 TMatrixD sideARMS(8,2);
1440 TMatrixD sideCRMS(8,2);
1441 TMatrixD sideACRMS(8,2);
1447 for (Int_t i=0; i<8;i++){
1448 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1449 his->Fit(&fg,"","", -0.15,0.15);
1450 sideARMS(i,0) = fg.GetParameter(2);
1451 his->SetDirectory(0);
1452 his->SetName(Form("Original SideA_%s",his->GetName()));
1453 his->SetTitle(Form("Original SideA_%s",his->GetTitle()));
1454 arrayOut.AddLast(his);
1455 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1456 his->Fit(&fg,"","", -0.15,0.15);
1457 sideARMS(i,1) = fg.GetParameter(2);
1458 his->SetDirectory(0);
1459 his->SetName(Form("Aligned SideA_%s",his->GetName()));
1460 his->SetTitle(Form("Aligned SideA_%s",his->GetTitle()));
1461 arrayOut.AddLast(his);
1467 for (Int_t i=0; i<8;i++){
1468 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1469 his->Fit(&fg,"","", -0.15,0.15);
1470 sideCRMS(i,0) = fg.GetParameter(2);
1471 his->SetDirectory(0);
1472 his->SetName(Form("Original SideC_%s",his->GetName()));
1473 his->SetTitle(Form("Original SideC_%s",his->GetTitle()));
1474 arrayOut.AddLast(his);
1475 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1476 his->Fit(&fg,"","", -0.15,0.15);
1477 sideCRMS(i,1) = fg.GetParameter(2);
1478 his->SetDirectory(0);
1479 his->SetName(Form("Aligned SideC_%s",his->GetName()));
1480 his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1481 arrayOut.AddLast(his);
1487 for (Int_t i=0; i<8;i++){
1488 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1489 his->Fit(&fg,"","", -0.15,0.15);
1490 sideACRMS(i,0) = fg.GetParameter(2);
1491 his->SetDirectory(0);
1492 his->SetName(Form("Original SideAC_%s",his->GetName()));
1493 his->SetTitle(Form("Original SideAC_%s",his->GetTitle()));
1494 arrayOut.AddLast(his);
1495 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1496 his->Fit(&fg,"","", -0.15,0.15);
1497 sideACRMS(i,1) = fg.GetParameter(2);
1498 his->SetDirectory(0);
1499 his->SetName(Form("Aligned SideC_%s",his->GetName()));
1500 his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1501 arrayOut.AddLast(his);
1509 void MakeFolders(TObjArray * arrayOut){
1513 TFolder *folderBase = new TFolder("TPC align","TPC align");
1516 TString maskDelta[8];
1517 TString maskAlign[2]={"Orig","Alig");
1518 for (Int_t i=0; i<8;i++){
1519 maskDelta[i] = kalmanFitNew->fLinearTrackDelta[i]->GetName();
1522 Int_t entries = arrayOut->GetEntries();
1529 void MergeKalman(const char * list = "kalmanFit.list"){
1535 TString currentFile;
1543 printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1544 if (currentFile.Length()==0) continue;
1545 TFile * ffit = TFile::Open(currentFile.Data());
1546 TEntryList * tlist = (TEntryList*) ffit->Get("eventList");
1547 TMatrixD * cMatrix = (TMatrixD*) ffit->Get("cutMarix");
1548 if (tlist&&cMatrix){
1549 printf("Track entries=%d\n",tlist->GetN());
1550 if (cMatrix->Sum()<=0.00000001) {
1551 printf("Problem with track selection\n");
1555 printf("Problem with track selection\n");
1560 AliTPCkalmanFit * fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFit");
1562 if (!kalmanFitNew) {kalmanFitNew= fit; continue;};
1563 kalmanFitNew->Add(fit);
1564 printf("Selected entries=%f\n",fit->fLinearTrackDelta[0]->GetEntries());
1572 currentFile.ReplaceAll("TPC","TPCOrig");
1573 printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1574 if (currentFile.Length()==0) continue;
1575 ffit = TFile::Open(currentFile.Data());
1576 fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFitOrig");
1578 if (!kalmanFitOrig) {kalmanFitOrig= fit; continue;};
1579 kalmanFitOrig->Add(fit);
1585 // save merged results
1587 TFile f("mergeKalmanFit.root","recreate");
1588 kalmanFitNew->Write("kalmanFitNew");
1589 kalmanFitOrig->Write("kalmanFitOrig");
1596 Example - submit alignment as batch jobs
1601 nclustersSkipOffset=0
1608 while [ $ntracksSkipOffset -lt $ntracksSkip ] ; do
1609 nclustersSkipOffset=0;
1610 while [ $nclustersSkipOffset -lt $nclustersSkip ] ; do
1611 echo Tr $ntracksSkipOffset Cl $nclustersSkipOffset;
1612 ver=kalmanDirTrack$ntracksSkipOffset$nclustersSkipOffset
1613 echo New Directory $ver
1616 cp $bDir/align.txt . ;
1617 ln -sf $bDir/kalmanFitApply.root . ;
1618 echo aliroot -q -b "$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C($ntracks,10000,0,$ndump,$ntracksSkip,$ntracksSkipOffset, $nclustersSkip,$nclustersSkipOffset,$isTest)" ;
1619 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)" ;
1620 nclustersSkipOffset=$(( $nclustersSkipOffset + 1 ))
1626 ntracksSkipOffset=$(( $ntracksSkipOffset + 1 ))
1627 echo Tr $ntracksSkipOffset;