AliRoot Core  ee782a0 (ee782a0)
CalibAlignKalman.C
Go to the documentation of this file.
1 
27 
28 #ifdef __CINT__
29 #else
30 #include <fstream>
31 #include "TSystem.h"
32 #include "TROOT.h"
33 #include "TRandom.h"
34 #include "TMath.h"
35 #include "TBits.h"
36 #include "TFormula.h"
37 #include "TF1.h"
38 #include "TLinearFitter.h"
39 #include "TFile.h"
40 #include "TChain.h"
41 #include "TCut.h"
42 #include "TEntryList.h"
43 #include "TH1F.h"
44 #include "THnSparse.h"
45 
46 
47 #include "AliSysInfo.h"
48 #include "AliExternalTrackParam.h"
49 #include "TTreeStream.h"
50 #include "AliTrackPointArray.h"
51 #include "AliLog.h"
52 #include "AliTPCTransformation.h"
53 #include "AliTPCkalmanFit.h"
54 #include "AliMathBase.h"
55 #include "AliXRDPROOFtoolkit.h"
56 #endif
57 
58 
59 //
60 //
61 // track point cuts
62 //
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
71 //
72 // mult. scattering cuts
73 /*
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
81 */
82 TMatrixD cutMatrix(4*7,2);
83 const Double_t rmsCut09[6]={0.15,0.9,0.015, 0.2, 0.9, 0.015};
84 
85 
86 
87 //
88 //
89 Int_t toSkip = 2; //
90 Int_t toSkipOffset = 0;
91 Int_t toSkipTrack = 2;
93 Int_t isFilterTest = 0;
94 //
95 // Track cuts
96 //
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};
103 
104 TCut cutAll="";
105 
106 
107 //
108 //
109 //
110 TChain *chainPoints=0;
111 TEntryList *elist=0;
115 
116 AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump);
117 void FilterTracks();
119 //
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);
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);
129 
130 
132 
134 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset);
135 
136 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump);
137 
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){
140 
142  toSkip=nSkip;
143  toSkipOffset= nSkipOffset;
144  toSkipTrack = nSkipTrack;
145  toSkipTrackOffset = nSkipTrackOffset;
146  isFilterTest = bfilterTest;
147  //
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();
155  }else{
156  printf("Not trnasformation specified\n");
157  }
158  //
159  //
160  //
161  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
162  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
163  //
164  AliXRDPROOFtoolkit tool;
165  chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0, maxFiles, startFile);
166  CalibAlignKalmanFit(npoints, trackDump);
167 }
168 
169 
170 
171 
172 AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump){
174 
176  FilterTracks();
177  kalmanFitNew = SetupFit();
178  kalmanFitOrig = SetupFit();
179  kalmanFitNew->Init();
180  kalmanFitOrig->Init();
181  return FitPointsLinear(maxTracks,trackDump);
182 }
183 
184 
185 
188 
189  AliTPCkalmanFit *kalmanFit = new AliTPCkalmanFit;
190  AddFitFieldCage(kalmanFit);
191  AddPhiScaling(kalmanFit);
192  AddDrift(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);
199  kalmanFit->Init();
200  return kalmanFit;
201 }
202 
203 
206 
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");
211  //
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);
216 
217  //
218 
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++){
222  if (!cA[ic]){
223  cA[ic]=new TCut;
224  *cA[ic]= *cSide[ic];
225  }
226  //
227  // cutP3
228  //
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);
234  //
235  // cutSP3
236  //
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);
242  //
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);
248  //
249  // cutP4
250  //
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);
261 
262  //
263  // cutM4
264  //
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);
270 
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);
276  //
277  //
278  //
279  cA[ic]=new TCut;
280  *cA[ic]= *cSide[ic]+*cP3[ic]+*cSP3[ic]+*cP4[ic]+*cM4[ic];
281  }
282  }
283  cutMatrix.Print();
284 
285  cutAll = (*cA[0])||(*cA[1])||(*cA[2])||(*cA[3])+"abs(mag)<0.01&&ncont>0&&p.fNPoints>120";
286 
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);
291 
292 
293 
294  chainPoints->Draw(">>listEL",cutAll,"entryList");
295  elist = (TEntryList*)gDirectory->Get("listEL");
296  chainPoints->SetEntryList(elist);
297  elist->SetDirectory(0);
298 }
299 
300 
301 
302 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump){
303  // create debug streamers
304 
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");
310  //
311  //
312  AliTrackPointArray *pointsNS=0;
313  Float_t mag=0;
314  Int_t time=0;
315  AliExternalTrackParam *param0=0;
316  AliExternalTrackParam *param1=0;
317  chainPoints->SetBranchAddress("p.",&pointsNS);
318  chainPoints->SetBranchAddress("p0In.",&param0);
319  chainPoints->SetBranchAddress("p1In.",&param1);
320  chainPoints->SetBranchAddress("mag",&mag);
321  chainPoints->SetBranchAddress("time",&time);
322  Int_t accepted=0;
323  printf("\n*\n*\n*Selected entries = %d\n*\n*\n*",Int_t(elist->GetN()));
324 
325  //
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;
331  //
332  AliTrackPointArray *points = AliTPCkalmanFit::SortPoints(*pointsNS);
333  if (kalmanFitApply) kalmanFitApply->ApplyCalibration(points,-1.);
334  //
335  // estimate and filter scattering
336  //
337  TVectorD *vecRMS09 = EstimateScatteringKalmanLinear(*points,*param0,*param1,pcstream);
338  if (!vecRMS09) continue;
339  Bool_t isOK=kTRUE;
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) {
347  delete points;
348  continue;
349  }
350  kalmanFitNew->PropagateTime(time);
351  //
352  //
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++;
361  //
362  if (accepted%50==0) {
363  kalmanFitNew->FitTrackLinear(*pointsF, pcstream);
364  }else{
365  if (idir==0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
366  if (idir!=0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
367  }
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);
372  }
373  AliSysInfo::AddStamp("trackFit", accepted,itrack);
374  delete pointsF;
375  delete spointsF;
376  }
377  delete points;
378  }
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");
385 
386  delete pcstream;
387  delete pcstreamOrig;
388  return kalmanFitNew;
389 }
390 
391 void QAPointsLinear(Int_t maxTracks, Int_t trackDump){
394 
395  // create debug streeamers
396  TTreeSRedirector *pcstreamNonCalib = new TTreeSRedirector("kalmanfitTPCQANonCalib.root");
397  TTreeSRedirector *pcstreamCalib = new TTreeSRedirector("kalmanfitTPCQACalib.root");
399  AliTPCkalmanFit *kalmanFitters[6]={0,0,0,0,0,0};
400  for (Int_t i=0;i<6;i++){
401  kalmanFitters[i]=SetupFit();
402  }
403  //
404  AliTrackPointArray *points=0;
405  AliExternalTrackParam *param0=0;
406  AliExternalTrackParam *param1=0;
407  Float_t mag=0;
408  Int_t time=0;
409  chainPoints->SetBranchAddress("p.",&points);
410  chainPoints->SetBranchAddress("mag",&mag);
411  chainPoints->SetBranchAddress("time",&time);
412  chainPoints->SetBranchAddress("p0In.",&param0);
413  chainPoints->SetBranchAddress("p1In.",&param1);
414 
415  Int_t accepted=0;
416  //
417 
418 
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;
424  //
425 
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;
434  //
435  if (idir==0) accepted++;
436  kalmanFitters[iscalib*3+idir+1]->DumpTrackLinear(*pointsF,0);
437  EstimateScatteringKalmanLinear(*pointsF,*param0,*param1,pcstream);
438  delete pointsF;
439  }
440  }
441  if (accepted%trackDump==0) {
442  printf("%d\n", accepted);
443  }
444  }
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;
454 }
455 
456 
457 void TestScattering(Int_t maxTracks, Int_t trackDump){
462 
463  TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPCMS.root");
464  //
465  //
466  AliTrackPointArray *points=0;
467  AliExternalTrackParam *param0=0;
468  AliExternalTrackParam *param1=0;
469  Float_t mag=0;
470  Int_t time=0;
471  chainPoints->SetBranchAddress("p.",&points);
472  chainPoints->SetBranchAddress("mag",&mag);
473  chainPoints->SetBranchAddress("time",&time);
474  chainPoints->SetBranchAddress("p0In.",&param0);
475  chainPoints->SetBranchAddress("p1In.",&param1);
476  Int_t accepted=0;
477  //
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;
482  //
483  AliTrackPointArray *pointsSorted = AliTPCkalmanFit::SortPoints(*points);
484  EstimateScatteringKalmanLinear(*pointsSorted,*param0,*param1,pcstream);
485  accepted++;
486  if (accepted%trackDump==0) {
487  printf("%d\n", accepted);
488  }
489  delete pointsSorted;
490  }
491  delete pcstream;
492 }
493 
494 
495 
496 
497 
498 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset){
500 
501  Int_t npoints = points.GetNPoints();
502  Int_t npointsF = (npoints-nskipOffset-1)/nskip;
503  AliTrackPoint point;
504  AliTrackPointArray *pointsF= new AliTrackPointArray(npointsF);
505  Int_t used=0;
506  for (Int_t ipoint=nskipOffset; ipoint<npoints; ipoint+=nskip){
507  //
508  if (!points.GetPoint(point,ipoint)) continue;
509  pointsF->AddPoint(used,&point);
510  used++;
511  if (used==npointsF) break;
512  }
513  return pointsF;
514 }
515 
516 
517 
520 
521  TLinearFitter lfitY(2,"pol1");
522  TLinearFitter lfitZ(2,"pol1");
523  TVectorD vecZ(2);
524  TVectorD vecY(2);
525  //
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);
533  //
534  // 1.b Fit the track in the rotated frame - MakeSeed
535  //
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);
547  }
548  if (TMath::Abs(maxX-minX)<kArmCut) return 0;
549  if (lfitY.GetNpoints()<knclCut) return 0;
550  //
551  lfitY.Eval();
552  lfitZ.Eval();
553  lfitY.GetParameters(vecY);
554  lfitZ.GetParameters(vecZ);
555  //
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;
560  //
561  //
562  Int_t accepted=0;
563  Int_t toBeUsed =0;
564  AliTrackPoint point;
565  AliTrackPointArray *pointsF=0;
566  for (Int_t iter=0; iter<2;iter++){
567  for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
568  //
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;
581  accepted++;
582  if (pointsF) pointsF->AddPoint(toBeUsed,&point);
583  toBeUsed++;
584  }
585  if (pcstream){
586  (*pcstream)<<"filter"<<
587  "iter="<<iter<<
588  "accepted="<<accepted<<
589  "minX="<<minX<<
590  "maxX="<<maxX<<
591  "vY.="<<&vecY<<
592  "vZ.="<<&vecZ<<
593  "chi2Y="<<chi2Y<<
594  "chi2Z="<<chi2Z<<
595  "\n";
596  }
597  if (accepted<knclCut) break;
598  if (iter==0) pointsF = new AliTrackPointArray(toBeUsed);
599  accepted=0;
600  toBeUsed=0;
601  }
602  return pointsF;
603 }
604 
605 
610 
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);
616  //
617  // 1. sort the points
618  //
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];
623  }
624  TMath::Sort(npoints, rxvector,indexes,kFALSE);
625  AliTrackPoint point;
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);
630  }
631  delete [] rxvector;
632  delete [] indexes;
633  return pointsSorted;
634 }
635 
639 
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.;
648 
649  Int_t npoints = points.GetNPoints();
650  if (npoints<kMinPoints) return 0;
651 
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);
655  }
656 
657  //
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);
661  //
662  //
663  TMatrixD trParamY(2,1),trCovarY(2,2);
664  TMatrixD trParamZ(2,1),trCovarZ(2,2);
665  Int_t kNmeas = 1;
666  Int_t nelem = 2;
667  //
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.;
682  matHk(0,0)=1;
683  matHk(0,1)=0;
684 
685  Double_t lastX = 0;
686  Int_t lastVolId=-1;
687  for (Int_t idir=0; idir<2;idir++){
688  // fit direction
689  //
690  for (Int_t ip=0; ip<npoints; ip++){
691  Int_t ipoint= 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];
697  //
698  if (ip==0){
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];
703  trParamY(0,0) = ry;
704  trParamY(1,0) = (rym-ry)/(rxm-rx);
705  trParamZ(0,0) = rz;
706  trParamZ(1,0) = (rzm-rz)/(rxm-rx);
707  //
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));
712  lastX = rx;
713  lastVolId = volId;
714  }
715  //
716  // Propagate
717  //
718  if ((volId%36)<18 && (lastVolId%36)>=18){
719  // A - C side cross
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);
724  lastVolId=volId;
725  }
726  if (volId!=lastVolId){
727  // volumeID change
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);
732  lastVolId=volId;
733  }
734  //
735  // Propagate
736  //
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;
744  trCovarY=covXk2;
745  covXk2=matFk*trCovarZ*matFkT;
746  trCovarZ=covXk2;
747 
748  // multiple scattering
749  trCovarY(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
750  trCovarZ(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
751  lastX=rx;
752  //
753  // Update
754  //
755  for (Int_t coord=0; coord<2;coord++){
756  TMatrixD* pvecXk = (coord==0)? &trParamY: &trParamZ;
757  TMatrixD* pcovXk = (coord==0)? &trCovarY: &trCovarZ;
758  //
759  TMatrixD& vecXk = *pvecXk;
760  TMatrixD& covXk = *pcovXk;
761  measR(0,0) = (coord==0) ? kResY:kResZ;
762  vecZk(0,0) = (coord==0) ? ry:rz;
763  //
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
767  matSk.Invert();
768  matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
769  //
770  covXk2= (mat1-(matKk*matHk));
771  covXk3 = covXk2*covXk;
772  vecXk += matKk*vecYk; // updated vector
773  covXk = covXk3;
774  }
775  //
776  // store parameters
777  //
778  (*vecPos[0])[ipoint]=rx;
779  (*vecPos[1])[ipoint]=ry;
780  (*vecPos[2])[ipoint]=rz;
781 
782  (*vecPos[4*idir+0+3])[ipoint]=trParamY(0,0);
783  (*vecPos[4*idir+1+3])[ipoint]=trParamY(1,0);
784  //
785  (*vecPos[4*idir+2+3])[ipoint]=trParamZ(0,0);
786  (*vecPos[4*idir+3+3])[ipoint]=trParamZ(1,0);
787  }
788  }
789  //
790  //
791  //
792  TVectorD vec(npoints);
793  TVectorD rms(6);
794  TVectorD rms09(6); // robust RMS - fraction 0.9
795  TVectorD mean09(6); // robust RMS - fraction 0.9
796  Double_t meanR,rmsR;
797  Int_t npoints09 = Int_t(npoints*0.9);
798  //
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);
802  rms09[0]=rmsR;
803  mean09[0]=meanR;
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);
807  rms09[1]=rmsR;
808  mean09[1]=meanR;
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);
812  rms09[2]=rmsR;
813  mean09[2]=meanR;
814  //
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);
818  rms09[3]=rmsR;
819  mean09[3]=meanR;
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);
823  rms09[4]=rmsR;
824  mean09[4]=meanR;
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);
828  rms09[5]=rmsR;
829  mean09[5]=meanR;
830 
831 
832  (*pcstream)<<"kf"<<
833  "p.="<<&points<<
834  "p0.="<<&p0<<
835  "p1.="<<&p1<<
836  "rms.="<<&rms<<
837  "rms09.="<<&rms09<<
838  "mean09.="<<&mean09<<
839  "py.="<<&trParamY<<
840  "pz.="<<&trParamZ<<
841  "cy.="<<&trCovarY<<
842  "cz.="<<&trCovarY<<
843  "\n";
844  for (Int_t i=0;i<11;i++){
845  delete vecPos[i];
846  }
847  return new TVectorD(rms09);
848 }
849 
850 
851 
852 
853 
854 
855 
856 
857 
860 
861  TVectorD fpar(10);
862  AliTPCTransformation * transformation=0;
863  char tname[100];
864  //
865  // linear R scaling and shift
866  //
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++){
870  fpar[0]=ipolR;
871  fpar[1]=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);
877  //
878  }
879  }
880  //
881  //
882  //Inner field cage
883  for (Int_t iside=0; iside<=1; iside++)
884  for (Int_t ipol=0; ipol<3; ipol++){
885  fpar[0]=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);
890  }
891  //
892  //
893  //Outer field cage
894  for (Int_t iside=0; iside<=1; iside++)
895  for (Int_t ipol=0; ipol<3; ipol++){
896  fpar[0]=ipol;
897  //Outer field cage
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);
902  }
903 }
904 
905 
911 
912  TBits maskInner(72);
913  TBits maskOuter(72);
914  for (Int_t i=0;i<72;i++){
915  if (i<36){
916  maskInner[i]=kTRUE;
917  }
918  if (i>=36){
919  maskOuter[i]=kTRUE;
920  }
921  }
922  TVectorD fpar(10);
923  AliTPCTransformation * transformation=0;
924  fpar[0]=1;
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);
931  //
932 }
933 
934 void AddDrift(AliTPCkalmanFit *kalmanFit){
936 
937  TVectorD fpar(10);
938  AliTPCTransformation * transformation=0;
939  fpar[0]=1;
940  transformation = new AliTPCTransformation("tTPCscalingZDrift", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDr", 0);
941  transformation->SetParams(0,4.,1.,&fpar);
942  kalmanFit->AddCalibration(transformation);
943  //
944  transformation = new AliTPCTransformation("tTPCscalingZDriftGy", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDrGy", 0);
945  transformation->SetParams(0,0.2,0.0,&fpar);
946  kalmanFit->AddCalibration(transformation);
947 }
948 
949 
950 
951 
952 void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
954 
955  TVectorD fpar(10);
956  fpar[0]=0; fpar[1]=0; fpar[2]=0;
957  char tname[1000];
958  AliTPCTransformation * transformation=0;
959  //
960  //
961  //
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
966  //
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);
971  //
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);
976  //
977  //
978  //
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);
984  //
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);
989  }
990 
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
995  //
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);
1000  //
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);
1005  //
1006  //
1007  //
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);
1013  //
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);
1018  }
1019 
1020 }
1021 
1022 
1023 
1024 
1025 void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1027 
1028  TVectorD fpar(10);
1029  fpar[0]=0; fpar[1]=0; fpar[2]=0;
1030  char tname[1000];
1031  AliTPCTransformation * transformation=0;
1032  //
1033  //
1034  //
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
1039  //
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);
1044  //
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);
1049  //
1050  //
1051  //
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);
1057  //
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);
1062  }
1063 
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
1068  //
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);
1073  //
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);
1078  //
1079  //
1080  //
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);
1086  //
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);
1091  }
1092 }
1093 
1094 
1095 
1098 
1099  TVectorD fpar(10);
1100  AliTPCTransformation * transformation=0;
1101  TBits maskInnerA(72);
1102  TBits maskInnerC(72);
1103  for (Int_t i=0;i<72;i++){
1104  if (i<36){
1105  if (i%36<18) maskInnerA[i]=kTRUE;
1106  if (i%36>=18) maskInnerC[i]=kTRUE;
1107  }
1108  }
1109  //
1110  //
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);
1117  //
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);
1124 }
1125 
1128 
1129  TVectorD fpar(10);
1130  AliTPCTransformation * transformation=0;
1131  Int_t fixSector =4;
1132  //
1133  for (Int_t isec=0; isec<36;isec++){
1134  TBits mask(72);
1135  mask[isec]=kTRUE;
1136  char tname[1000];
1137  //
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);
1142  //
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);
1147 
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);
1152 
1153  }
1154  //
1155  for (Int_t isec=0; isec<36;isec++){
1156  if (isec%18==fixSector) continue;
1157  TBits mask(72);
1158  mask[isec] =kTRUE;
1159  mask[isec+36]=kTRUE;
1160  char tname[1000];
1161  //
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);
1166  //
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);
1171 
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);
1176  }
1177 
1178 
1179 }
1180 
1181 
1182 void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1184 
1185  TVectorD fpar(10);
1186  AliTPCTransformation * transformation=0;
1187 
1188  for (Int_t i=0; i<=ncos;i++){
1189  char tname[1000];
1190  fpar[0]=i; // cos frequency
1191  fpar[1]=0; // no sinus
1192  fpar[2]=1; // relative misalignment
1193  //
1194  // Local X shift
1195  //
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);
1204  //
1205  //
1206  // Local Y shift
1207  //
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);
1216  //
1217  //
1218  //
1219  // Z rotation
1220  //
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);
1229  }
1230  //
1231  for (Int_t i=1; i<=nsin;i++){
1232  char tname[1000];
1233  fpar[0]=0; // cos frequency
1234  fpar[1]=i; // sinus frequency
1235  fpar[2]=1; // relative misalignment
1236  //
1237  // Local X shift
1238  //
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);
1247  //
1248  //
1249  // Local Y shift
1250  //
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);
1259  //
1260  //
1261  //
1262  // Z rotation
1263  //
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);
1272  }
1273 }
1274 
1275 void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1277 
1278  TVectorD fpar(10);
1279  AliTPCTransformation * transformation=0;
1280 
1281  for (Int_t i=0; i<=ncos;i++){
1282  char tname[1000];
1283  fpar[0]=i; // cos frequency
1284  fpar[1]=0; // no sinus
1285  fpar[2]=0; // absolute misalignment
1286  if (i>0){
1287  //
1288  // A side is reference
1289  // local x
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);
1294  if (i>1){
1295  //
1296  // Local Y shift
1297  //
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);
1302  }
1303  //
1304  //
1305  }
1306  //
1307  // C side to vary
1308  // local x
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);
1313  //
1314  // Local Y shift
1315  //
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);
1320  //
1321  //
1322  // Z rotation - independent
1323  //
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);
1332  }
1333 
1334 
1335 
1336 
1337  //
1338  //
1339  //
1340  for (Int_t i=1; i<=nsin;i++){
1341  char tname[1000];
1342  fpar[0]=0; // non cos frequency
1343  fpar[1]=i; // sinus frequency
1344  fpar[2]=0; // absolute misalignment
1345  //
1346  // Local X shift
1347  //
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);
1356  //
1357  //
1358  // Local Y shift
1359  //
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);
1368  //
1369  //
1370  //
1371  // Z rotation
1372  //
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);
1381  }
1382 }
1383 
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);
1390  }
1391 }
1392 
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);
1399  }
1400 }
1401 
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);
1408  }
1409 }
1410 
1411 
1412 void DumpQA1D( TObjArray &arrayOut){
1414 
1415  TF1 fg("fg","gaus");
1416  TMatrixD sideARMS(8,2);
1417  TMatrixD sideCRMS(8,2);
1418  TMatrixD sideACRMS(8,2);
1419  TH1 *his=0;
1420  //
1421  // A side
1422  //
1423  SelectNonPixelA();
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);
1439  }
1440  //
1441  // C side
1442  //
1443  SelectNonPixelC();
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);
1459  }
1460  //
1461  // AC side
1462  //
1463  SelectPixel();
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);
1479  }
1480  printf("DumQA\n");
1481  sideARMS.Print();
1482  sideCRMS.Print();
1483  sideACRMS.Print();
1484 }
1485 
1486 void MakeFolders(TObjArray * arrayOut){
1488 
1489  TFolder *folderBase = new TFolder("TPC align","TPC align");
1490  //
1491  //
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();
1496  }
1497 
1498  Int_t entries = arrayOut->GetEntries();
1499  //
1500 }
1501 
1502 
1503 
1504 
1505 void MergeKalman(const char * list = "kalmanFit.list"){
1507 
1508  ifstream in;
1509  in.open(list);
1510  TString currentFile;
1511  kalmanFitNew= 0;
1512  Int_t counter=0;
1513  while(in.good()) {
1514  //
1515  // calibrated
1516  //
1517  in >> 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");
1527  continue;
1528  }
1529  }else{
1530  printf("Problem with track selection\n");
1531  continue;
1532  }
1533  //
1534 
1535  AliTPCkalmanFit * fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFit");
1536  if (!fit) continue;
1537  if (!kalmanFitNew) {kalmanFitNew= fit; continue;};
1538  kalmanFitNew->Add(fit);
1539  printf("Selected entries=%f\n",fit->fLinearTrackDelta[0]->GetEntries());
1540  //delete tlist;
1541  //delete cMatrix;
1542  delete fit;
1543  delete ffit;
1544  //
1545  // original
1546  //
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");
1552  if (!fit) continue;
1553  if (!kalmanFitOrig) {kalmanFitOrig= fit; continue;};
1554  kalmanFitOrig->Add(fit);
1555  delete fit;
1556  delete ffit;
1557  counter++;
1558  }
1559  //
1560  // save merged results
1561  //
1562  TFile f("mergeKalmanFit.root","recreate");
1563  kalmanFitNew->Write("kalmanFitNew");
1564  kalmanFitOrig->Write("kalmanFitOrig");
1565  f.Close();
1566 }
1567 
1568 
1569 
1570 /*
1571  Example - submit alignment as batch jobs
1572  ifile=0;
1573  ntracksSkip=200
1574  ntracksSkipOffset=0
1575  nclustersSkip=2
1576  nclustersSkipOffset=0
1577  ntracks=1000000
1578  ndump=5
1579  isTest=0
1580  bDir=`pwd`
1581  ver=aaa;
1582 
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
1589  mkdir $ver;
1590  cd $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 ))
1596  cd $bDir;
1597  echo $bDir;
1598  done;
1599  cd $bDir;
1600  echo $bDir;
1601  ntracksSkipOffset=$(( $ntracksSkipOffset + 1 ))
1602  echo Tr $ntracksSkipOffset;
1603 done
1604 
1605 
1606 
1607 */
1608 
1609 
1610 
1611 
1612 
1613 
1614 
1615 
1616 
1617 
1618 
TCut cutAll
Bool_t AddPoint(Int_t i, const AliTrackPoint *p)
THnSparse * fLinearTrackDelta[12]
linear tracks matching residuals - delta
void AddPhiScaling(AliTPCkalmanFit *kalmanFit)
Int_t toSkipOffset
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TFile * Open(const char *filename, Long64_t &nevents)
TOOLKIT for chain manipulation:
const Float_t krmsZcutGlobal
void FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug=0, Float_t scalingRMSY=1., Float_t scalingRMSZ=1.)
TCut * cM4[4]
#define TObjArray
TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1, TTreeSRedirector *pcstream)
TVectorD vec
Definition: AnalyzeLaser.C:8
AliTPCkalmanFit * kalmanFitApply
void AddFitFieldCage(AliTPCkalmanFit *kalmanFit)
AliTrackPointArray * SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset)
TFile * GetFile()
Definition: TTreeStream.h:95
void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin)
void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin)
const Double_t kNsigma
const Float_t * GetY() const
Int_t GetNPoints() const
void DumpQA1D(TObjArray &arrayOut)
void TestScattering(Int_t maxTracks, Int_t trackDump)
TROOT * gROOT
void SelectNonPixelA()
AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump)
const Float_t krmsYcutGlobal
TTreeSRedirector * pcstream
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)
TCut * cA[4]
void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin)
TMatrixD cutMatrix(4 *7, 2)
void SelectNonPixelC()
const Float_t krmsYcut
npoints
Definition: driftITSTPC.C:85
AliSplineFit fit
Definition: CalibTime.C:30
void PropagateTime(Int_t time)
const UShort_t * GetVolumeID() const
TCut * cSP3[4]
void ApplyCalibration(AliTrackPointArray *array, Double_t csign)
Int_t toSkipTrackOffset
TCut * cP4[4]
void SelectPixel()
Bool_t GetPoint(AliTrackPoint &p, Int_t i) const
TEntryList * elist
Int_t toSkip
const Double_t kArmCut
void SetParams(Double_t param, Double_t sigma, Double_t sigma2Time, const TVectorD *const fixedParams)
const Float_t krmsZcut
AliTrackPointArray * FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream)
Should represent general non linear transformation.
TCut * cP3[4]
const Float_t kSigmaCut
void FilterTracks()
const Int_t knclCut
AliTPCkalmanFit * kalmanFitNew
const Double_t rmsCut09[6]
Bool_t DumpCalib(const char *mask=0, Float_t correlationCut=-1)
TF1 * f
Definition: interpolTest.C:21
AliTPCkalmanFit * kalmanFitOrig
void DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug)
void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit)
Int_t toSkipTrack
AliTrackPointArray * SortPoints(AliTrackPointArray &points)
void QAPointsLinear(Int_t maxTracks, Int_t trackDump)
static AliTrackPointArray * SortPoints(AliTrackPointArray &points)
void AddCalibration(AliTPCTransformation *calib)
static Int_t BuildBasicFormulas()
static AliMagF::BMap_t mag
void AddDrift(AliTPCkalmanFit *kalmanFit)
TCut * cSide[4]
void MakeFolders(TObjArray *arrayOut)
static void AddStamp(const char *sname, Int_t id0=-1, Int_t id1=-1, Int_t id2=-1, Int_t id3=-1)
Definition: AliSysInfo.cxx:179
void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin)
class TVectorT< Double_t > TVectorD
static TBits * BitsSide(Bool_t aside)
void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit)
TChain * chainPoints
AliTPCkalmanFit * SetupFit()
static TChain * MakeChainRandom(const char *fileIn, const char *treeName, const char *fName=0, Int_t maxFiles=-1, Int_t startFile=0, Int_t checkLevel=0)
static void EvaluateUni(const Int_t nvectors, const Double_t *data, Double_t &mean, Double_t &sigma, const Int_t hSub)
Definition: AliMathBase.cxx:60
Int_t isFilterTest
void Add(const AliTPCkalmanFit *kalman)
AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump)
const Float_t * GetZ() const
class TMatrixT< Double_t > TMatrixD
const Float_t * GetX() const