]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/CalibAlignKalman.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibAlignKalman.C
CommitLineData
35d81915 1/// \file CalibAlignKalman.C
2///
3/// ~~~{.cpp}
4/// gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
5/// gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
6/// gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
7///
8/// gROOT->LoadMacro("$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C+");
9///
10/// AliTPCTransformation::BuildBasicFormulas();
11///
12/// AliXRDPROOFtoolkit tool;
13/// chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0,50000);
14/// chainPoints->Lookup();
15///
16/// chainMS = tool.MakeChainRandom("kalmanFit.list","kf",0,50000);
17/// chainMS->Lookup();
18///
19/// chainFP = tool.MakeChainRandom("kalmanFit.list","filter",0,50000);
20/// chainFP->Lookup();
21///
22/// CalibAlignKalmanFit(40000,1);
23/// kalmanFit0->DumpCorelation(0.8);
24/// TFile f("kalmanfitTPC.root");
25/// ~~~
171ba885 26
27
d8b91e72 28#ifdef __CINT__
29#else
171ba885 30#include <fstream>
171ba885 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"
d8b91e72 44#include "THnSparse.h"
171ba885 45
46
d8b91e72 47#include "AliSysInfo.h"
48#include "AliExternalTrackParam.h"
171ba885 49#include "TTreeStream.h"
50#include "AliTrackPointArray.h"
51#include "AliLog.h"
52#include "AliTPCTransformation.h"
53#include "AliTPCkalmanFit.h"
d8b91e72 54#include "AliMathBase.h"
171ba885 55#include "AliXRDPROOFtoolkit.h"
d8b91e72 56#endif
57
58
59//
60//
61// track point cuts
62//
63const Float_t krmsYcut = 0.2; // cluster RMS cut in Y - residuals form local fit
64const Float_t krmsZcut = 0.2; // cluster RMS cut in Z - residuals from local fit
65const Float_t krmsYcutGlobal= 0.2; // cluster RMS cut in Y - residuals form global fit
66const Float_t krmsZcutGlobal= 2.0; // cluster RMS cut in Z - residuals from global fit
67const Float_t kSigmaCut = 5.; // clusters to be removed
68const Int_t knclCut = 80; // minimal number of clusters
69const Double_t kArmCut = 50.; // minimal level arm
70const Double_t kNsigma = 5.; // cut on track level
71//
72// mult. scattering cuts
73/*
74TCut cutClY("rms09.fElements[0]<0.15");
75TCut cuttY("rms09.fElements[1]/rms09.fElements[0]<0.9");
76TCut cutkY("rms09.fElements[2]<0.015");
77TCut cutClZ("rms09.fElements[3]<0.2");
78TCut cuttZ("rms09.fElements[4]/rms09.fElements[3]<0.9");
79TCut cutkZ("rms09.fElements[5]<0.015");
80TCut cutMS=cutClY+cuttY+cutkY+cutClZ+cuttZ+cutkZ
81*/
82TMatrixD cutMatrix(4*7,2);
83const Double_t rmsCut09[6]={0.15,0.9,0.015, 0.2, 0.9, 0.015};
84
85
86
87//
88//
89Int_t toSkip = 2; //
90Int_t toSkipOffset = 0;
91Int_t toSkipTrack = 2;
92Int_t toSkipTrackOffset = 0;
93Int_t isFilterTest = 0;
94//
95// Track cuts
96//
97TCut * cSide[4]={0,0,0,0};
98TCut *cP3[4]={0,0,0,0};
99TCut *cSP3[4]={0,0,0,0};
100TCut *cP4[4]={0,0,0,0};
101TCut *cM4[4]={0,0,0,0};
102TCut *cA[4]={0,0,0,0};
103
104TCut cutAll="";
171ba885 105
106
d8b91e72 107//
108//
171ba885 109//
110TChain *chainPoints=0;
111TEntryList *elist=0;
d8b91e72 112AliTPCkalmanFit * kalmanFitNew=0;
113AliTPCkalmanFit * kalmanFitOrig=0;
114AliTPCkalmanFit * kalmanFitApply=0;
171ba885 115
d8b91e72 116AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump);
171ba885 117void FilterTracks();
118AliTPCkalmanFit * SetupFit();
d8b91e72 119//
171ba885 120void AddFitFieldCage(AliTPCkalmanFit *kalmanFit);
121void AddPhiScaling(AliTPCkalmanFit *kalmanFit);
d8b91e72 122void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin );
123void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
171ba885 124void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit);
d8b91e72 125void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit);
126void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
127void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
128void AddDrift(AliTPCkalmanFit *kalmanFit);
129
130
171ba885 131AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream);
171ba885 132
d8b91e72 133TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream);
134AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset);
135
136AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump);
137
138void 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){
35d81915 139 ///
140
d8b91e72 141 AliTPCTransformation::BuildBasicFormulas();
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 }
171ba885 158 //
159 //
160 //
161 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
162 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
d8b91e72 163 //
171ba885 164 AliXRDPROOFtoolkit tool;
165 chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0, maxFiles, startFile);
d8b91e72 166 CalibAlignKalmanFit(npoints, trackDump);
171ba885 167}
168
169
170
171ba885 171
d8b91e72 172AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump){
35d81915 173 /// Fitting procedure
174
171ba885 175 AliTPCTransformation::BuildBasicFormulas();
176 FilterTracks();
d8b91e72 177 kalmanFitNew = SetupFit();
178 kalmanFitOrig = SetupFit();
179 kalmanFitNew->Init();
180 kalmanFitOrig->Init();
181 return FitPointsLinear(maxTracks,trackDump);
171ba885 182}
183
184
185
186AliTPCkalmanFit * SetupFit(){
35d81915 187 ///
188
171ba885 189 AliTPCkalmanFit *kalmanFit = new AliTPCkalmanFit;
190 AddFitFieldCage(kalmanFit);
191 AddPhiScaling(kalmanFit);
d8b91e72 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();
171ba885 200 return kalmanFit;
201}
202
203
204void FilterTracks(){
35d81915 205 ///
d8b91e72 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
171ba885 294 chainPoints->Draw(">>listEL",cutAll,"entryList");
295 elist = (TEntryList*)gDirectory->Get("listEL");
296 chainPoints->SetEntryList(elist);
297 elist->SetDirectory(0);
298}
299
300
301
d8b91e72 302AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump){
35d81915 303 // create debug streamers
304
171ba885 305 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPC.root");
d8b91e72 306 TTreeSRedirector *pcstreamOrig = new TTreeSRedirector("kalmanfitTPCOrig.root");
307 pcstream->GetFile()->cd();
308 cutMatrix.Write("cutMarix");
309 elist->Write("eventList");
171ba885 310 //
311 //
d8b91e72 312 AliTrackPointArray *pointsNS=0;
171ba885 313 Float_t mag=0;
314 Int_t time=0;
d8b91e72 315 AliExternalTrackParam *param0=0;
316 AliExternalTrackParam *param1=0;
317 chainPoints->SetBranchAddress("p.",&pointsNS);
318 chainPoints->SetBranchAddress("p0In.",&param0);
319 chainPoints->SetBranchAddress("p1In.",&param1);
171ba885 320 chainPoints->SetBranchAddress("mag",&mag);
321 chainPoints->SetBranchAddress("time",&time);
322 Int_t accepted=0;
d8b91e72 323 printf("\n*\n*\n*Selected entries = %d\n*\n*\n*",Int_t(elist->GetN()));
324
171ba885 325 //
326 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
d8b91e72 327 if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
171ba885 328 Int_t entry=chainPoints->GetEntryNumber(itrack);
329 chainPoints->GetEntry(entry);
171ba885 330 if (accepted>maxTracks) break;
d8b91e72 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 //
171ba885 353 for (Int_t idir=-1; idir<=1; idir++){
d8b91e72 354 AliTrackPointArray *pointsF = FilterPoints(*points,idir, pcstream);
171ba885 355 if (!pointsF) continue;
d8b91e72 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);
171ba885 360 if (idir==0) accepted++;
361 //
362 if (accepted%50==0) {
d8b91e72 363 kalmanFitNew->FitTrackLinear(*pointsF, pcstream);
171ba885 364 }else{
d8b91e72 365 if (idir==0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
366 if (idir!=0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
171ba885 367 }
d8b91e72 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);
171ba885 374 delete pointsF;
d8b91e72 375 delete spointsF;
171ba885 376 }
d8b91e72 377 delete points;
171ba885 378 }
379 pcstream->GetFile()->cd();
d8b91e72 380 kalmanFitNew->Write("kalmanFit");
381 pcstreamOrig->GetFile()->cd();
382 kalmanFitOrig->Write("kalmanFitOrig");
383 pcstreamOrig->GetFile()->cd();
384 if (kalmanFitApply) kalmanFitApply->Write("kalmanFitApply");
385
171ba885 386 delete pcstream;
d8b91e72 387 delete pcstreamOrig;
388 return kalmanFitNew;
171ba885 389}
390
d8b91e72 391void QAPointsLinear(Int_t maxTracks, Int_t trackDump){
35d81915 392 /// check the consistency of kalman fit
393 /// Apply transformation
394
d8b91e72 395 // create debug streeamers
396 TTreeSRedirector *pcstreamNonCalib = new TTreeSRedirector("kalmanfitTPCQANonCalib.root");
397 TTreeSRedirector *pcstreamCalib = new TTreeSRedirector("kalmanfitTPCQACalib.root");
398 TTreeSRedirector *pcstream=0;
399 AliTPCkalmanFit *kalmanFitters[6]={0,0,0,0,0,0};
400 for (Int_t i=0;i<6;i++){
401 kalmanFitters[i]=SetupFit();
402 }
403 //
d8b91e72 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 //
171ba885 417
418
d8b91e72 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 //
171ba885 425
d8b91e72 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
457void TestScattering(Int_t maxTracks, Int_t trackDump){
35d81915 458 /// test Multiple scattering algorithm
459 /// Apply transformation
460 ///
461 /// create debug streeamers
462
d8b91e72 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;
171ba885 477 //
d8b91e72 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
498AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset){
35d81915 499 /// create new array with skipped points
500
d8b91e72 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
518AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream){
35d81915 519 /// Filter points - input points for KalmanFilter
520
171ba885 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();
d8b91e72 558 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
559 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
171ba885 560 //
561 //
562 Int_t accepted=0;
d8b91e72 563 Int_t toBeUsed =0;
171ba885 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;
171ba885 581 accepted++;
d8b91e72 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";
171ba885 596 }
597 if (accepted<knclCut) break;
d8b91e72 598 if (iter==0) pointsF = new AliTrackPointArray(toBeUsed);
171ba885 599 accepted=0;
d8b91e72 600 toBeUsed=0;
171ba885 601 }
602 return pointsF;
603}
604
605
d8b91e72 606AliTrackPointArray * SortPoints(AliTrackPointArray &points){
35d81915 607 /// Creates the array - points sorted according radius - neccessay for kalman fit
608 ///
609 /// 0. choose the frame - rotation angle
610
d8b91e72 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
636TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream){
35d81915 637 /// Algorithm - 0. Fit the track forward and backward
638 /// - 1. Store the current parameters in each point
d8b91e72 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
171ba885 856
857
858void AddFitFieldCage(AliTPCkalmanFit *kalmanFit){
35d81915 859 /// Add radial scaling due field cage
860
171ba885 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
906void AddPhiScaling(AliTPCkalmanFit *kalmanFit){
35d81915 907 /// Add linear local phi scaling -
908 /// separate IROC/OROC - A side/C side
909 /// "tscalingLocalPhiIROC"
910 /// "tscalingLocalPhiOROC"
d8b91e72 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
934void AddDrift(AliTPCkalmanFit *kalmanFit){
35d81915 935 /// Add drift velocity transformation
936
171ba885 937 TVectorD fpar(10);
938 AliTPCTransformation * transformation=0;
939 fpar[0]=1;
d8b91e72 940 transformation = new AliTPCTransformation("tTPCscalingZDrift", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDr", 0);
941 transformation->SetParams(0,4.,1.,&fpar);
171ba885 942 kalmanFit->AddCalibration(transformation);
943 //
d8b91e72 944 transformation = new AliTPCTransformation("tTPCscalingZDriftGy", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDrGy", 0);
945 transformation->SetParams(0,0.2,0.0,&fpar);
946 kalmanFit->AddCalibration(transformation);
171ba885 947}
948
d8b91e72 949
950
951
952void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
35d81915 953 ///
954
171ba885 955 TVectorD fpar(10);
d8b91e72 956 fpar[0]=0; fpar[1]=0; fpar[2]=0;
957 char tname[1000];
171ba885 958 AliTPCTransformation * transformation=0;
171ba885 959 //
171ba885 960 //
171ba885 961 //
d8b91e72 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
171ba885 1020}
1021
d8b91e72 1022
1023
1024
1025void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
35d81915 1026 /// z tilting absolute (sector) and relative (IROC-OROC)
1027
171ba885 1028 TVectorD fpar(10);
d8b91e72 1029 fpar[0]=0; fpar[1]=0; fpar[2]=0;
1030 char tname[1000];
171ba885 1031 AliTPCTransformation * transformation=0;
171ba885 1032 //
171ba885 1033 //
1034 //
d8b91e72 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 }
171ba885 1092}
1093
1094
1095
1096void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit){
35d81915 1097 ///
1098
171ba885 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
d8b91e72 1126void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit){
35d81915 1127 ///
1128
d8b91e72 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
1182void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
35d81915 1183 ///
1184
d8b91e72 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
1275void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
35d81915 1276 ///
1277
d8b91e72 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
1384void SelectPixel(){
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
1393void SelectNonPixelA(){
1394 for (Int_t i=0;i<8;i++){
1395 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1396 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1397 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1398 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1399 }
1400}
1401
1402void SelectNonPixelC(){
1403 for (Int_t i=0;i<8;i++){
1404 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1405 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1406 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1407 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1408 }
1409}
1410
1411
1412void DumpQA1D( TObjArray &arrayOut){
35d81915 1413 ///
1414
d8b91e72 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
1486void MakeFolders(TObjArray * arrayOut){
35d81915 1487 ///
1488
d8b91e72 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
171ba885 1505void MergeKalman(const char * list = "kalmanFit.list"){
35d81915 1506 ///
1507
171ba885 1508 ifstream in;
1509 in.open(list);
1510 TString currentFile;
d8b91e72 1511 kalmanFitNew= 0;
171ba885 1512 Int_t counter=0;
1513 while(in.good()) {
d8b91e72 1514 //
1515 // calibrated
1516 //
1517 in >> currentFile;
171ba885 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());
d8b91e72 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
171ba885 1535 AliTPCkalmanFit * fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFit");
1536 if (!fit) continue;
d8b91e72 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);
171ba885 1555 delete fit;
1556 delete ffit;
1557 counter++;
1558 }
d8b91e72 1559 //
1560 // save merged results
1561 //
1562 TFile f("mergeKalmanFit.root","recreate");
1563 kalmanFitNew->Write("kalmanFitNew");
1564 kalmanFitOrig->Write("kalmanFitOrig");
1565 f.Close();
171ba885 1566}
1567
1568
d8b91e72 1569
171ba885 1570/*
d8b91e72 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
171ba885 1580 bDir=`pwd`
d8b91e72 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;
1603done
171ba885 1604
171ba885 1605
1606
d8b91e72 1607*/
1608
171ba885 1609
1610
1611
1612
1613
1614
1615
1616
1617
1618