]> 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
Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibAlignKalman.C
CommitLineData
171ba885 1/*
2
d8b91e72 3 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
4 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
5 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
171ba885 6
d8b91e72 7 gROOT->LoadMacro("$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C+");
171ba885 8
d8b91e72 9 AliTPCTransformation::BuildBasicFormulas();
171ba885 10
d8b91e72 11 AliXRDPROOFtoolkit tool;
12 chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0,50000);
13 chainPoints->Lookup();
14 //
15 chainMS = tool.MakeChainRandom("kalmanFit.list","kf",0,50000);
16 chainMS->Lookup();
17
18 chainFP = tool.MakeChainRandom("kalmanFit.list","filter",0,50000);
19 chainFP->Lookup();
20
21// CalibAlignKalmanFit(40000,1);
171ba885 22// kalmanFit0->DumpCorelation(0.8);
23// TFile f("kalmanfitTPC.root");
d8b91e72 24*/
171ba885 25
26
d8b91e72 27#ifdef __CINT__
28#else
171ba885 29#include <fstream>
171ba885 30#include "TSystem.h"
31#include "TROOT.h"
32#include "TRandom.h"
33#include "TMath.h"
34#include "TBits.h"
35#include "TFormula.h"
36#include "TF1.h"
37#include "TLinearFitter.h"
38#include "TFile.h"
39#include "TChain.h"
40#include "TCut.h"
41#include "TEntryList.h"
42#include "TH1F.h"
d8b91e72 43#include "THnSparse.h"
171ba885 44
45
d8b91e72 46#include "AliSysInfo.h"
47#include "AliExternalTrackParam.h"
171ba885 48#include "TTreeStream.h"
49#include "AliTrackPointArray.h"
50#include "AliLog.h"
51#include "AliTPCTransformation.h"
52#include "AliTPCkalmanFit.h"
d8b91e72 53#include "AliMathBase.h"
171ba885 54#include "AliXRDPROOFtoolkit.h"
d8b91e72 55#endif
56
57
58//
59//
60// track point cuts
61//
62const Float_t krmsYcut = 0.2; // cluster RMS cut in Y - residuals form local fit
63const Float_t krmsZcut = 0.2; // cluster RMS cut in Z - residuals from local fit
64const Float_t krmsYcutGlobal= 0.2; // cluster RMS cut in Y - residuals form global fit
65const Float_t krmsZcutGlobal= 2.0; // cluster RMS cut in Z - residuals from global fit
66const Float_t kSigmaCut = 5.; // clusters to be removed
67const Int_t knclCut = 80; // minimal number of clusters
68const Double_t kArmCut = 50.; // minimal level arm
69const Double_t kNsigma = 5.; // cut on track level
70//
71// mult. scattering cuts
72/*
73TCut cutClY("rms09.fElements[0]<0.15");
74TCut cuttY("rms09.fElements[1]/rms09.fElements[0]<0.9");
75TCut cutkY("rms09.fElements[2]<0.015");
76TCut cutClZ("rms09.fElements[3]<0.2");
77TCut cuttZ("rms09.fElements[4]/rms09.fElements[3]<0.9");
78TCut cutkZ("rms09.fElements[5]<0.015");
79TCut cutMS=cutClY+cuttY+cutkY+cutClZ+cuttZ+cutkZ
80*/
81TMatrixD cutMatrix(4*7,2);
82const Double_t rmsCut09[6]={0.15,0.9,0.015, 0.2, 0.9, 0.015};
83
84
85
86//
87//
88Int_t toSkip = 2; //
89Int_t toSkipOffset = 0;
90Int_t toSkipTrack = 2;
91Int_t toSkipTrackOffset = 0;
92Int_t isFilterTest = 0;
93//
94// Track cuts
95//
96TCut * cSide[4]={0,0,0,0};
97TCut *cP3[4]={0,0,0,0};
98TCut *cSP3[4]={0,0,0,0};
99TCut *cP4[4]={0,0,0,0};
100TCut *cM4[4]={0,0,0,0};
101TCut *cA[4]={0,0,0,0};
102
103TCut cutAll="";
171ba885 104
105
d8b91e72 106//
107//
171ba885 108//
109TChain *chainPoints=0;
110TEntryList *elist=0;
d8b91e72 111AliTPCkalmanFit * kalmanFitNew=0;
112AliTPCkalmanFit * kalmanFitOrig=0;
113AliTPCkalmanFit * kalmanFitApply=0;
171ba885 114
d8b91e72 115AliTPCkalmanFit * CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump);
171ba885 116void FilterTracks();
117AliTPCkalmanFit * SetupFit();
d8b91e72 118//
171ba885 119void AddFitFieldCage(AliTPCkalmanFit *kalmanFit);
120void AddPhiScaling(AliTPCkalmanFit *kalmanFit);
d8b91e72 121void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin );
122void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
171ba885 123void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit);
d8b91e72 124void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit);
125void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
126void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
127void AddDrift(AliTPCkalmanFit *kalmanFit);
128
129
171ba885 130AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream);
171ba885 131
d8b91e72 132TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream);
133AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset);
134
135AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump);
136
137void 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){
138 //
139 //
140 //
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){
171ba885 173 //
d8b91e72 174 // Fitting procedure
171ba885 175 //
176 AliTPCTransformation::BuildBasicFormulas();
177 FilterTracks();
d8b91e72 178 kalmanFitNew = SetupFit();
179 kalmanFitOrig = SetupFit();
180 kalmanFitNew->Init();
181 kalmanFitOrig->Init();
182 return FitPointsLinear(maxTracks,trackDump);
171ba885 183}
184
185
186
187AliTPCkalmanFit * SetupFit(){
188 //
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(){
205 //
206 //
207 //
d8b91e72 208
209 cSide[0] = new TCut("cutAA","p0In.fP[1]>0&&p1In.fP[1]>0");
210 cSide[1] = new TCut("cutCC","p0In.fP[1]<0&&p1In.fP[1]<0");
211 cSide[2] = new TCut("cutAC","p0In.fP[1]>0&&p1In.fP[1]<0");
212 cSide[3] = new TCut("cutCA","p0In.fP[1]<0&&p1In.fP[1]>0");
213 //
214 TH1F * phisP3 = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
215 TH1F * phisSP3 = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
216 TH1F * phisP4 = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
217 TH1F * phisM4 = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
218
219 //
220
221 TF1 *fg = new TF1("fg","gaus");
222 for (Int_t iter=0; iter<2; iter++){
223 for (Int_t ic=0;ic<4;ic++){
224 if (!cA[ic]){
225 cA[ic]=new TCut;
226 *cA[ic]= *cSide[ic];
227 }
228 //
229 // cutP3
230 //
231 chainPoints->Draw("p0In.fP[3]+p1In.fP[3]>>hhisP3",*cA[ic],"goff");
232 phisP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisP3->GetMaximumBin())-0.003,phisP3->GetBinCenter(phisP3->GetMaximumBin())+0.003);
233 cP3[ic]=new TCut(Form("abs(p0In.fP[3]+p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
234 cutMatrix(7*ic+0,0) = fg->GetParameter(1);
235 cutMatrix(7*ic+0,1) = fg->GetParameter(2);
236 //
237 // cutSP3
238 //
239 chainPoints->Draw("p0Out.fP[3]-p0In.fP[3]>>hhisSP3",*cA[ic],"goff");
240 phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
241 cSP3[ic]=new TCut(Form("abs(p0Out.fP[3]-p0In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
242 cutMatrix(7*ic+1,0) = fg->GetParameter(1);
243 cutMatrix(7*ic+1,1) = fg->GetParameter(2);
244 //
245 chainPoints->Draw("p1Out.fP[3]-p1In.fP[3]>>hhisSP3",*cA[ic],"goff");
246 phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
247 *cSP3[ic]+=Form("abs(p1Out.fP[3]-p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
248 cutMatrix(7*ic+2,0) = fg->GetParameter(1);
249 cutMatrix(7*ic+2,1) = fg->GetParameter(2);
250 //
251 // cutP4
252 //
253 chainPoints->Draw("p0Out.fP[4]>>hhisP4",*cA[ic],"goff");
254 phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
255 cP4[ic]=new TCut(Form("abs(p0Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
256 cutMatrix(7*ic+3,0) = fg->GetParameter(1);
257 cutMatrix(7*ic+3,1) = fg->GetParameter(2);
258 chainPoints->Draw("p1Out.fP[4]>>hhisP4",*cA[ic],"goff");
259 phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
260 *cP4[ic]+=Form("abs(p1Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
261 cutMatrix(7*ic+4,0) = fg->GetParameter(1);
262 cutMatrix(7*ic+4,1) = fg->GetParameter(2);
263
264 //
265 // cutM4
266 //
267 chainPoints->Draw("p0Out.fP[4]-p0In.fP[4]>>hhisM4",*cA[ic],"goff");
268 phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
269 cM4[ic]=new TCut(Form("abs(p0Out.fP[4]-p0In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
270 cutMatrix(7*ic+5,0) = fg->GetParameter(1);
271 cutMatrix(7*ic+5,1) = fg->GetParameter(2);
272
273 chainPoints->Draw("p1Out.fP[4]-p1In.fP[4]>>hhisM4",*cA[ic],"goff");
274 phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
275 *cM4[ic]+=Form("abs(p1Out.fP[4]-p1In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
276 cutMatrix(7*ic+6,0) = fg->GetParameter(1);
277 cutMatrix(7*ic+6,1) = fg->GetParameter(2);
278 //
279 //
280 //
281 cA[ic]=new TCut;
282 *cA[ic]= *cSide[ic]+*cP3[ic]+*cSP3[ic]+*cP4[ic]+*cM4[ic];
283 }
284 }
285 cutMatrix.Print();
286
287 cutAll = (*cA[0])||(*cA[1])||(*cA[2])||(*cA[3])+"abs(mag)<0.01&&ncont>0&&p.fNPoints>120";
288
289 delete phisP3; // = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
290 delete phisSP3; // = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
291 delete phisP4;// = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
292 delete phisM4;// = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
293
294
295
171ba885 296 chainPoints->Draw(">>listEL",cutAll,"entryList");
297 elist = (TEntryList*)gDirectory->Get("listEL");
298 chainPoints->SetEntryList(elist);
299 elist->SetDirectory(0);
300}
301
302
303
d8b91e72 304AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump){
171ba885 305 //
306 //
307 //
308 // create debug streeamers
309 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPC.root");
d8b91e72 310 TTreeSRedirector *pcstreamOrig = new TTreeSRedirector("kalmanfitTPCOrig.root");
311 pcstream->GetFile()->cd();
312 cutMatrix.Write("cutMarix");
313 elist->Write("eventList");
171ba885 314 //
315 //
d8b91e72 316 AliTrackPointArray *pointsNS=0;
171ba885 317 Float_t mag=0;
318 Int_t time=0;
d8b91e72 319 AliExternalTrackParam *param0=0;
320 AliExternalTrackParam *param1=0;
321 chainPoints->SetBranchAddress("p.",&pointsNS);
322 chainPoints->SetBranchAddress("p0In.",&param0);
323 chainPoints->SetBranchAddress("p1In.",&param1);
171ba885 324 chainPoints->SetBranchAddress("mag",&mag);
325 chainPoints->SetBranchAddress("time",&time);
326 Int_t accepted=0;
d8b91e72 327 printf("\n*\n*\n*Selected entries = %d\n*\n*\n*",Int_t(elist->GetN()));
328
171ba885 329 //
330 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
d8b91e72 331 if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
171ba885 332 Int_t entry=chainPoints->GetEntryNumber(itrack);
333 chainPoints->GetEntry(entry);
171ba885 334 if (accepted>maxTracks) break;
d8b91e72 335 //
336 AliTrackPointArray *points = AliTPCkalmanFit::SortPoints(*pointsNS);
337 if (kalmanFitApply) kalmanFitApply->ApplyCalibration(points,-1.);
338 //
339 // estimate and filter scattering
340 //
341 TVectorD *vecRMS09 = EstimateScatteringKalmanLinear(*points,*param0,*param1,pcstream);
342 if (!vecRMS09) continue;
343 Bool_t isOK=kTRUE;
344 if ((*vecRMS09)[0] >rmsCut09[0]) isOK=kFALSE;
345 if ((*vecRMS09)[1]/(*vecRMS09)[0] >rmsCut09[1]) isOK=kFALSE;;
346 if ((*vecRMS09)[2] >rmsCut09[2]) isOK=kFALSE;
347 if ((*vecRMS09)[3] >rmsCut09[3]) isOK=kFALSE;
348 if ((*vecRMS09)[4]/(*vecRMS09)[0] >rmsCut09[4]) isOK=kFALSE;
349 if ((*vecRMS09)[5] >rmsCut09[5]) isOK=kFALSE;
350 if (!isOK || isFilterTest) {
351 delete points;
352 continue;
353 }
354 kalmanFitNew->PropagateTime(time);
355 //
356 //
171ba885 357 for (Int_t idir=-1; idir<=1; idir++){
d8b91e72 358 AliTrackPointArray *pointsF = FilterPoints(*points,idir, pcstream);
171ba885 359 if (!pointsF) continue;
d8b91e72 360 AliTrackPointArray *spointsF = 0;
361 // we skip points for alignemnt but not for QA
362 if (idir==0) spointsF = SkipPoints(*pointsF, toSkip*2, toSkipOffset);
363 if (idir!=0) spointsF = SkipPoints(*pointsF, toSkip, toSkipOffset);
171ba885 364 if (idir==0) accepted++;
365 //
366 if (accepted%50==0) {
d8b91e72 367 kalmanFitNew->FitTrackLinear(*pointsF, pcstream);
171ba885 368 }else{
d8b91e72 369 if (idir==0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
370 if (idir!=0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
171ba885 371 }
d8b91e72 372 if (idir==0) kalmanFitNew->DumpTrackLinear(*pointsF,pcstream);
373 if (idir==0) kalmanFitOrig->DumpTrackLinear(*pointsF,pcstreamOrig);
374 if (accepted%trackDump==0) {
375 printf("%d\n", accepted);
376 }
377 AliSysInfo::AddStamp("trackFit", accepted,itrack);
171ba885 378 delete pointsF;
d8b91e72 379 delete spointsF;
171ba885 380 }
d8b91e72 381 delete points;
171ba885 382 }
383 pcstream->GetFile()->cd();
d8b91e72 384 kalmanFitNew->Write("kalmanFit");
385 pcstreamOrig->GetFile()->cd();
386 kalmanFitOrig->Write("kalmanFitOrig");
387 pcstreamOrig->GetFile()->cd();
388 if (kalmanFitApply) kalmanFitApply->Write("kalmanFitApply");
389
171ba885 390 delete pcstream;
d8b91e72 391 delete pcstreamOrig;
392 return kalmanFitNew;
171ba885 393}
394
d8b91e72 395void QAPointsLinear(Int_t maxTracks, Int_t trackDump){
396 //
397 // check the consistency of kalman fit
398 // Apply transformation
399 //
400 // create debug streeamers
401 TTreeSRedirector *pcstreamNonCalib = new TTreeSRedirector("kalmanfitTPCQANonCalib.root");
402 TTreeSRedirector *pcstreamCalib = new TTreeSRedirector("kalmanfitTPCQACalib.root");
403 TTreeSRedirector *pcstream=0;
404 AliTPCkalmanFit *kalmanFitters[6]={0,0,0,0,0,0};
405 for (Int_t i=0;i<6;i++){
406 kalmanFitters[i]=SetupFit();
407 }
408 //
409 //
410 AliTrackPointArray *points=0;
411 AliExternalTrackParam *param0=0;
412 AliExternalTrackParam *param1=0;
413 Float_t mag=0;
414 Int_t time=0;
415 chainPoints->SetBranchAddress("p.",&points);
416 chainPoints->SetBranchAddress("mag",&mag);
417 chainPoints->SetBranchAddress("time",&time);
418 chainPoints->SetBranchAddress("p0In.",&param0);
419 chainPoints->SetBranchAddress("p1In.",&param1);
420
421 Int_t accepted=0;
422 //
171ba885 423
424
d8b91e72 425 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
426 if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
427 Int_t entry=chainPoints->GetEntryNumber(itrack);
428 chainPoints->GetEntry(entry);
429 if (accepted>maxTracks) break;
430 //
171ba885 431
d8b91e72 432 AliTrackPointArray pointsCalib(*points);
433 for (Int_t iscalib=0; iscalib<1;iscalib++){
434 if (iscalib>0) kalmanFitNew->ApplyCalibration(&pointsCalib,-1.);
435 if (iscalib==0) pcstream=pcstreamNonCalib;
436 if (iscalib>0) pcstream=pcstreamCalib;
437 for (Int_t idir=-1; idir<=1; idir++){
438 AliTrackPointArray *pointsF = FilterPoints(pointsCalib,idir, pcstream);
439 if (!pointsF) continue;
440 //
441 if (idir==0) accepted++;
442 kalmanFitters[iscalib*3+idir+1]->DumpTrackLinear(*pointsF,0);
443 EstimateScatteringKalmanLinear(*pointsF,*param0,*param1,pcstream);
444 delete pointsF;
445 }
446 }
447 if (accepted%trackDump==0) {
448 printf("%d\n", accepted);
449 }
450 }
451 pcstreamCalib->GetFile()->cd();
452 kalmanFitters[0]->Write("fitUpNonCalib");
453 kalmanFitters[1]->Write("fitUpDownNonCalib");
454 kalmanFitters[2]->Write("fitDownNonCalib");
455 kalmanFitters[3]->Write("fitUpCalib");
456 kalmanFitters[4]->Write("fitUpDownCalib");
457 kalmanFitters[5]->Write("fitDownCalib");
458 delete pcstreamCalib;
459 delete pcstreamNonCalib;
460}
461
462
463void TestScattering(Int_t maxTracks, Int_t trackDump){
464 //
465 // test Multiple scattering algorithm
466 // Apply transformation
467 //
468 // create debug streeamers
469 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTPCMS.root");
470 //
471 //
472 AliTrackPointArray *points=0;
473 AliExternalTrackParam *param0=0;
474 AliExternalTrackParam *param1=0;
475 Float_t mag=0;
476 Int_t time=0;
477 chainPoints->SetBranchAddress("p.",&points);
478 chainPoints->SetBranchAddress("mag",&mag);
479 chainPoints->SetBranchAddress("time",&time);
480 chainPoints->SetBranchAddress("p0In.",&param0);
481 chainPoints->SetBranchAddress("p1In.",&param1);
482 Int_t accepted=0;
171ba885 483 //
d8b91e72 484 for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
485 Int_t entry=chainPoints->GetEntryNumber(itrack);
486 chainPoints->GetEntry(entry);
487 if (accepted>maxTracks) break;
488 //
489 AliTrackPointArray *pointsSorted = AliTPCkalmanFit::SortPoints(*points);
490 EstimateScatteringKalmanLinear(*pointsSorted,*param0,*param1,pcstream);
491 accepted++;
492 if (accepted%trackDump==0) {
493 printf("%d\n", accepted);
494 }
495 delete pointsSorted;
496 }
497 delete pcstream;
498}
499
500
501
502
503
504AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset){
171ba885 505 //
d8b91e72 506 // create new array with skipped points
507 //
508 Int_t npoints = points.GetNPoints();
509 Int_t npointsF = (npoints-nskipOffset-1)/nskip;
510 AliTrackPoint point;
511 AliTrackPointArray *pointsF= new AliTrackPointArray(npointsF);
512 Int_t used=0;
513 for (Int_t ipoint=nskipOffset; ipoint<npoints; ipoint+=nskip){
514 //
515 if (!points.GetPoint(point,ipoint)) continue;
516 pointsF->AddPoint(used,&point);
517 used++;
518 if (used==npointsF) break;
519 }
520 return pointsF;
521}
522
523
524
525AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream){
526 //
527 // Filter points - input points for KalmanFilter
528 //
171ba885 529 //
530 TLinearFitter lfitY(2,"pol1");
531 TLinearFitter lfitZ(2,"pol1");
532 TVectorD vecZ(2);
533 TVectorD vecY(2);
534 //
535 lfitY.StoreData(kTRUE);
536 lfitZ.StoreData(kTRUE);
537 Int_t npoints = points.GetNPoints();
538 if (npoints<2) return 0;
539 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
540 Double_t ca = TMath::Cos(currentAlpha);
541 Double_t sa = TMath::Sin(currentAlpha);
542 //
543 // 1.b Fit the track in the rotated frame - MakeSeed
544 //
545 Double_t maxX =-10000, minX=10000;
546 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
547 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
548 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
549 Double_t rz = points.GetZ()[ipoint];
550 if (dir== 1 && rx<0) continue;
551 if (dir==-1 && rx>0) continue;
552 if (maxX<rx) maxX=rx;
553 if (minX>rx) minX=rx;
554 lfitY.AddPoint(&rx,ry,1);
555 lfitZ.AddPoint(&rx,rz,1);
556 }
557 if (TMath::Abs(maxX-minX)<kArmCut) return 0;
558 if (lfitY.GetNpoints()<knclCut) return 0;
559 //
560 lfitY.Eval();
561 lfitZ.Eval();
562 lfitY.GetParameters(vecY);
563 lfitZ.GetParameters(vecZ);
564 //
565 Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
566 Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
d8b91e72 567 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
568 if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
171ba885 569 //
570 //
571 Int_t accepted=0;
d8b91e72 572 Int_t toBeUsed =0;
171ba885 573 AliTrackPoint point;
574 AliTrackPointArray *pointsF=0;
575 for (Int_t iter=0; iter<2;iter++){
576 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
577 //
578 if (!points.GetPoint(point,ipoint)) continue;
579 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
580 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
581 Double_t rz = points.GetZ()[ipoint];
582 if (dir== 1 && rx<0) continue;
583 if (dir==-1 && rx>0) continue;
584 Double_t erry = TMath::Sqrt(chi2Y);
585 Double_t errz = TMath::Sqrt(chi2Z);
586 Double_t fy = vecY[0]+vecY[1]*rx;
587 Double_t fz = vecZ[0]+vecZ[1]*rx;
588 if (TMath::Abs(fy-ry)>erry*kSigmaCut) continue;
589 if (TMath::Abs(fz-rz)>errz*kSigmaCut) continue;
171ba885 590 accepted++;
d8b91e72 591 if (pointsF) pointsF->AddPoint(toBeUsed,&point);
592 toBeUsed++;
593 }
594 if (pcstream){
595 (*pcstream)<<"filter"<<
596 "iter="<<iter<<
597 "accepted="<<accepted<<
598 "minX="<<minX<<
599 "maxX="<<maxX<<
600 "vY.="<<&vecY<<
601 "vZ.="<<&vecZ<<
602 "chi2Y="<<chi2Y<<
603 "chi2Z="<<chi2Z<<
604 "\n";
171ba885 605 }
606 if (accepted<knclCut) break;
d8b91e72 607 if (iter==0) pointsF = new AliTrackPointArray(toBeUsed);
171ba885 608 accepted=0;
d8b91e72 609 toBeUsed=0;
171ba885 610 }
611 return pointsF;
612}
613
614
d8b91e72 615AliTrackPointArray * SortPoints(AliTrackPointArray &points){
616 //
617 //Creates the array - points sorted according radius - neccessay for kalman fit
618 //
619 //
620 // 0. choose the frame - rotation angle
621 //
622 Int_t npoints = points.GetNPoints();
623 if (npoints<1) return 0;
624 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
625 Double_t ca = TMath::Cos(currentAlpha);
626 Double_t sa = TMath::Sin(currentAlpha);
627 //
628 // 1. sort the points
629 //
630 Double_t *rxvector = new Double_t[npoints];
631 Int_t *indexes = new Int_t[npoints];
632 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
633 rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
634 }
635 TMath::Sort(npoints, rxvector,indexes,kFALSE);
636 AliTrackPoint point;
637 AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
638 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
639 if (!points.GetPoint(point,indexes[ipoint])) continue;
640 pointsSorted->AddPoint(ipoint,&point);
641 }
642 delete [] rxvector;
643 delete [] indexes;
644 return pointsSorted;
645}
646
647TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream){
648 //
649 // Algorithm - 0. Fit the track forward and backward
650 // - 1. Store the current parameters in each point
651
652 const Int_t kMinPoints= 70;
653 const Double_t kResY = 0.1;
654 const Double_t kResZ = 0.1;
655 const Double_t kMisY = 0.02;
656 const Double_t kMisZ = 0.02;
657 const Double_t kLArm = 120.;
658 const Double_t kACsideFac = 10.;
659 const Double_t kMS0 = kResY/20.;
660
661 Int_t npoints = points.GetNPoints();
662 if (npoints<kMinPoints) return 0;
663
664 TVectorD *vecPos[11]={0,0,0,0,0,0,0,0,0,0,0};
665 for (Int_t i=0;i<11;i++){
666 vecPos[i]=new TVectorD(npoints);
667 }
668
669 //
670 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
671 Double_t ca = TMath::Cos(currentAlpha);
672 Double_t sa = TMath::Sin(currentAlpha);
673 //
674 //
675 TMatrixD trParamY(2,1),trCovarY(2,2);
676 TMatrixD trParamZ(2,1),trCovarZ(2,2);
677 Int_t kNmeas = 1;
678 Int_t nelem = 2;
679 //
680 TMatrixD matHk(kNmeas,nelem); // vector to mesurement
681 TMatrixD matHkT(nelem,kNmeas); // helper matrix Hk transpose
682 TMatrixD matFk(nelem,nelem);
683 TMatrixD matFkT(nelem,nelem);
684 TMatrixD vecYk(kNmeas,1); // Innovation or measurement residual
685 TMatrixD vecZk(kNmeas,1); // Innovation or measurement residual
686 TMatrixD measR(kNmeas,kNmeas);
687 TMatrixD matSk(kNmeas,kNmeas); // Innovation (or residual) covariance
688 TMatrixD matKk(nelem,kNmeas); // Optimal Kalman gain
689 TMatrixD covXk2(nelem,nelem); // helper matrix
690 TMatrixD covXk3(nelem,nelem); // helper matrix
691 TMatrixD mat1(nelem,nelem);
692 mat1(0,0)=1.; mat1(0,1)=0.;
693 mat1(1,0)=0.; mat1(1,1)=1.;
694 matHk(0,0)=1;
695 matHk(0,1)=0;
696
697 Double_t lastX = 0;
698 Int_t lastVolId=-1;
699 for (Int_t idir=0; idir<2;idir++){
700 // fit direction
701 //
702 for (Int_t ip=0; ip<npoints; ip++){
703 Int_t ipoint= ip;
704 if (idir>0) ipoint = npoints-ip-1;
705 Double_t rx = ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
706 Double_t ry = -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
707 Double_t rz = points.GetZ()[ipoint];
708 Int_t volId= points.GetVolumeID()[ipoint];
709 //
710 if (ip==0){
711 // set initital parameters and covariance - use first and middle point
712 Double_t rxm = ca*points.GetX()[npoints/2]+sa*points.GetY()[npoints/2];
713 Double_t rym = -sa*points.GetX()[npoints/2]+ca*points.GetY()[npoints/2];
714 Double_t rzm = points.GetZ()[npoints/2];
715 trParamY(0,0) = ry;
716 trParamY(1,0) = (rym-ry)/(rxm-rx);
717 trParamZ(0,0) = rz;
718 trParamZ(1,0) = (rzm-rz)/(rxm-rx);
719 //
720 trCovarY(0,0) = kResY*kResY;
721 trCovarY(1,1) = (kResY*kResY)/((rxm-rx)*(rxm-rx));
722 trCovarZ(0,0) = kResZ*kResZ;
723 trCovarZ(1,1) = (kResZ*kResZ)/((rxm-rx)*(rxm-rx));
724 lastX = rx;
725 lastVolId = volId;
726 }
727 //
728 // Propagate
729 //
730 if ((volId%36)<18 && (lastVolId%36)>=18){
731 // A - C side cross
732 trCovarY(0,0)+=kMisY*kMisY*kACsideFac;
733 trCovarZ(0,0)+=kMisZ*kMisZ*kACsideFac;
734 trCovarY(1,1)+=kACsideFac*(kMisY*kMisY)/(kLArm*kLArm);;
735 trCovarZ(1,1)+=kACsideFac*(kMisZ*kMisZ)/(kLArm*kLArm);
736 lastVolId=volId;
737 }
738 if (volId!=lastVolId){
739 // volumeID change
740 trCovarY(0,0)+=kMisY*kMisY;
741 trCovarZ(0,0)+=kMisZ*kMisZ;
742 trCovarY(1,1)+=(kMisY*kMisY)/(kLArm*kLArm);;
743 trCovarZ(1,1)+=(kMisZ*kMisZ)/(kLArm*kLArm);
744 lastVolId=volId;
745 }
746 //
747 // Propagate
748 //
749 Double_t deltaX=rx-lastX;
750 trParamY(0,0)+=deltaX*trParamY(1,0);
751 trParamZ(0,0)+=deltaX*trParamZ(1,0);
752 matFk(0,0)=1; matFk(0,1)=deltaX;
753 matFk(1,0)=0; matFk(1,1)=1;
754 matFkT=matFk.T(); matFk.T();
755 covXk2=matFk*trCovarY*matFkT;
756 trCovarY=covXk2;
757 covXk2=matFk*trCovarZ*matFkT;
758 trCovarZ=covXk2;
759
760 // multiple scattering
761 trCovarY(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
762 trCovarZ(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
763 lastX=rx;
764 //
765 // Update
766 //
767 for (Int_t coord=0; coord<2;coord++){
768 TMatrixD* pvecXk = (coord==0)? &trParamY: &trParamZ;
769 TMatrixD* pcovXk = (coord==0)? &trCovarY: &trCovarZ;
770 //
771 TMatrixD& vecXk = *pvecXk;
772 TMatrixD& covXk = *pcovXk;
773 measR(0,0) = (coord==0) ? kResY:kResZ;
774 vecZk(0,0) = (coord==0) ? ry:rz;
775 //
776 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
777 matHkT=matHk.T(); matHk.T();
778 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
779 matSk.Invert();
780 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
781 //
782 covXk2= (mat1-(matKk*matHk));
783 covXk3 = covXk2*covXk;
784 vecXk += matKk*vecYk; // updated vector
785 covXk = covXk3;
786 }
787 //
788 // store parameters
789 //
790 (*vecPos[0])[ipoint]=rx;
791 (*vecPos[1])[ipoint]=ry;
792 (*vecPos[2])[ipoint]=rz;
793
794 (*vecPos[4*idir+0+3])[ipoint]=trParamY(0,0);
795 (*vecPos[4*idir+1+3])[ipoint]=trParamY(1,0);
796 //
797 (*vecPos[4*idir+2+3])[ipoint]=trParamZ(0,0);
798 (*vecPos[4*idir+3+3])[ipoint]=trParamZ(1,0);
799 }
800 }
801 //
802 //
803 //
804 TVectorD vec(npoints);
805 TVectorD rms(6);
806 TVectorD rms09(6); // robust RMS - fraction 0.9
807 TVectorD mean09(6); // robust RMS - fraction 0.9
808 Double_t meanR,rmsR;
809 Int_t npoints09 = Int_t(npoints*0.9);
810 //
811 vec=(*(vecPos[3])-*(vecPos[1]));
812 rms[0]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms cluster y
813 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
814 rms09[0]=rmsR;
815 mean09[0]=meanR;
816 vec=(*(vecPos[7])-*(vecPos[3]));
817 rms[1]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track y
818 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
819 rms09[1]=rmsR;
820 mean09[1]=meanR;
821 vec=(*(vecPos[8])-*(vecPos[4]));
822 rms[2]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track ky
823 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
824 rms09[2]=rmsR;
825 mean09[2]=meanR;
826 //
827 vec=(*(vecPos[5])-*(vecPos[2]));
828 rms[3]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms cluster z
829 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
830 rms09[3]=rmsR;
831 mean09[3]=meanR;
832 vec=(*(vecPos[9])-*(vecPos[5]));
833 rms[4]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track z
834 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
835 rms09[4]=rmsR;
836 mean09[4]=meanR;
837 vec=(*(vecPos[10])-*(vecPos[6]));
838 rms[5]=TMath::RMS(npoints, vec.GetMatrixArray()); // rms track kz
839 AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
840 rms09[5]=rmsR;
841 mean09[5]=meanR;
842
843
844 (*pcstream)<<"kf"<<
845 "p.="<<&points<<
846 "p0.="<<&p0<<
847 "p1.="<<&p1<<
848 "rms.="<<&rms<<
849 "rms09.="<<&rms09<<
850 "mean09.="<<&mean09<<
851 "py.="<<&trParamY<<
852 "pz.="<<&trParamZ<<
853 "cy.="<<&trCovarY<<
854 "cz.="<<&trCovarY<<
855 "\n";
856 for (Int_t i=0;i<11;i++){
857 delete vecPos[i];
858 }
859 return new TVectorD(rms09);
860}
861
862
863
864
865
866
867
171ba885 868
869
870void AddFitFieldCage(AliTPCkalmanFit *kalmanFit){
871 //
872 // Add radial scaling due field cage
873 //
874 TVectorD fpar(10);
875 AliTPCTransformation * transformation=0;
876 char tname[100];
877 //
878 // linear R scaling and shift
879 //
880 for (Int_t iside=0; iside<=1; iside++)
881 for (Int_t ipolR=0; ipolR<2; ipolR++){
882 for (Int_t ipolZ=0; ipolZ<3; ipolZ++){
883 fpar[0]=ipolR;
884 fpar[1]=ipolZ;
885 if (ipolR+ipolZ==0) continue;
886 sprintf(tname,"tTPCscalingRPolR%dDr%dSide%d",ipolR,ipolZ,iside);
887 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0, 1);
888 transformation->SetParams(0,0.2,0,&fpar);
889 kalmanFit->AddCalibration(transformation);
890 //
891 }
892 }
893 //
894 //
895 //Inner field cage
896 for (Int_t iside=0; iside<=1; iside++)
897 for (Int_t ipol=0; ipol<3; ipol++){
898 fpar[0]=ipol;
899 sprintf(tname,"tTPCscalingRIFC%dSide%d",ipol,iside);
900 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRIFC",0,0, 1);
901 transformation->SetParams(0,0.2,0,&fpar);
902 kalmanFit->AddCalibration(transformation);
903 }
904 //
905 //
906 //Outer field cage
907 for (Int_t iside=0; iside<=1; iside++)
908 for (Int_t ipol=0; ipol<3; ipol++){
909 fpar[0]=ipol;
910 //Outer field cage
911 sprintf(tname,"tTPCscalingROFC%dSide%d",ipol,iside);
912 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingROFC",0,0, 1);
913 transformation->SetParams(0,0.2,0,&fpar);
914 kalmanFit->AddCalibration(transformation);
915 }
916}
917
918
919void AddPhiScaling(AliTPCkalmanFit *kalmanFit){
920 //
d8b91e72 921 // Add linear local phi scaling -
922 // separate IROC/OROC - A side/C side
923 // "tscalingLocalPhiIROC"
924 // "tscalingLocalPhiOROC"
925
926 TBits maskInner(72);
927 TBits maskOuter(72);
928 for (Int_t i=0;i<72;i++){
929 if (i<36){
930 maskInner[i]=kTRUE;
931 }
932 if (i>=36){
933 maskOuter[i]=kTRUE;
934 }
935 }
936 TVectorD fpar(10);
937 AliTPCTransformation * transformation=0;
938 fpar[0]=1;
939 transformation = new AliTPCTransformation("tscalingLocalPhiIROC", new TBits(maskInner), 0,"TPCscalingPhiLocal",0, 1);
940 transformation->SetParams(0,0.02,0,&fpar);
941 kalmanFit->AddCalibration(transformation);
942 transformation = new AliTPCTransformation("tscalingLocalPhiOROC", new TBits(maskOuter), 0,"TPCscalingPhiLocal",0, 1);
943 transformation->SetParams(0,0.02,0,&fpar);
944 kalmanFit->AddCalibration(transformation);
945 //
946}
947
948void AddDrift(AliTPCkalmanFit *kalmanFit){
949 //
950 // Add drift velocity transformation
171ba885 951 //
952 TVectorD fpar(10);
953 AliTPCTransformation * transformation=0;
954 fpar[0]=1;
d8b91e72 955 transformation = new AliTPCTransformation("tTPCscalingZDrift", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDr", 0);
956 transformation->SetParams(0,4.,1.,&fpar);
171ba885 957 kalmanFit->AddCalibration(transformation);
958 //
d8b91e72 959 transformation = new AliTPCTransformation("tTPCscalingZDriftGy", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDrGy", 0);
960 transformation->SetParams(0,0.2,0.0,&fpar);
961 kalmanFit->AddCalibration(transformation);
171ba885 962}
963
d8b91e72 964
965
966
967void AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
968 //
969 //
970 //
171ba885 971 TVectorD fpar(10);
d8b91e72 972 fpar[0]=0; fpar[1]=0; fpar[2]=0;
973 char tname[1000];
171ba885 974 AliTPCTransformation * transformation=0;
171ba885 975 //
171ba885 976 //
171ba885 977 //
d8b91e72 978 for (Int_t i=0; i<=ncos;i++){
979 fpar[0]=i; // cos frequency
980 fpar[1]=0; // no sinus
981 fpar[2]=1; // relative misalignment
982 //
983 sprintf(tname,"tTPCDeltaZIROCOROCSideA_Cos%d",i);
984 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
985 transformation->SetParams(0,0.03,0,&fpar);
986 kalmanFit->AddCalibration(transformation);
987 //
988 sprintf(tname,"tTPCDeltaZIROCOROCSideC_Cos%d",i);
989 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
990 transformation->SetParams(0,0.03,0,&fpar);
991 kalmanFit->AddCalibration(transformation);
992 //
993 //
994 //
995 fpar[2]=0; // absolute misalignment
996 sprintf(tname,"tTPCDeltaZSectorSideA_Cos%d",i);
997 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
998 transformation->SetParams(0,0.1,0,&fpar);
999 kalmanFit->AddCalibration(transformation);
1000 //
1001 sprintf(tname,"tTPCDeltaZSectorSideC_Cos%d",i);
1002 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1003 transformation->SetParams(0,0.1,0,&fpar);
1004 kalmanFit->AddCalibration(transformation);
1005 }
1006
1007 for (Int_t i=1; i<=nsin;i++){
1008 fpar[0]=0; // cos frequency
1009 fpar[1]=i; // sinus frequncy
1010 fpar[2]=1; // relative misalignment
1011 //
1012 sprintf(tname,"tTPCDeltaZIROCOROCSideA_Sin%d",i);
1013 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
1014 transformation->SetParams(0,0.03,0,&fpar);
1015 kalmanFit->AddCalibration(transformation);
1016 //
1017 sprintf(tname,"tTPCDeltaZIROCOROCSideC_Sin%d",i);
1018 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1019 transformation->SetParams(0,0.03,0,&fpar);
1020 kalmanFit->AddCalibration(transformation);
1021 //
1022 //
1023 //
1024 fpar[2]=0; // absolute misalignment
1025 sprintf(tname,"tTPCDeltaZSectorSideA_Sin%d",i);
1026 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" , 0);
1027 transformation->SetParams(0,0.1,0,&fpar);
1028 kalmanFit->AddCalibration(transformation);
1029 //
1030 sprintf(tname,"tTPCDeltaZSectorSideC_Sin%d",i);
1031 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" , 0);
1032 transformation->SetParams(0,0.1,0,&fpar);
1033 kalmanFit->AddCalibration(transformation);
1034 }
1035
171ba885 1036}
1037
d8b91e72 1038
1039
1040
1041void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
171ba885 1042 //
d8b91e72 1043 // z tilting absolute (sector) and relative (IROC-OROC)
171ba885 1044 //
1045 //
1046 TVectorD fpar(10);
d8b91e72 1047 fpar[0]=0; fpar[1]=0; fpar[2]=0;
1048 char tname[1000];
171ba885 1049 AliTPCTransformation * transformation=0;
171ba885 1050 //
171ba885 1051 //
1052 //
d8b91e72 1053 for (Int_t i=0; i<=ncos;i++){
1054 fpar[0]=i; // cos frequency
1055 fpar[1]=0; // sinus frequency
1056 fpar[2]=1; // relative misalignment
1057 //
1058 sprintf(tname,"tTPCTiltingZIROCOROCSideA_Cos%d",i);
1059 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1060 transformation->SetParams(0,0.03,0,&fpar);
1061 kalmanFit->AddCalibration(transformation);
1062 //
1063 sprintf(tname,"tTPCTiltingZIROCOROCSideC_Cos%d",i);
1064 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1065 transformation->SetParams(0,0.03,0,&fpar);
1066 kalmanFit->AddCalibration(transformation);
1067 //
1068 //
1069 //
1070 fpar[2]=0; // absolute misalignment
1071 sprintf(tname,"tTPCTiltingZSectorSideA_Cos%d",i);
1072 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1073 transformation->SetParams(0,0.1,0,&fpar);
1074 kalmanFit->AddCalibration(transformation);
1075 //
1076 sprintf(tname,"tTPCTiltingZSectorSideC_Cos%d",i);
1077 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1078 transformation->SetParams(0,0.1,0,&fpar);
1079 kalmanFit->AddCalibration(transformation);
1080 }
1081
1082 for (Int_t i=1; i<=nsin;i++){
1083 fpar[0]=0; // cos frequency
1084 fpar[1]=i; // sinus frequncy
1085 fpar[2]=1; // relative misalignment
1086 //
1087 sprintf(tname,"tTPCTiltingZIROCOROCSideA_Sin%d",i);
1088 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1089 transformation->SetParams(0,0.03,0,&fpar);
1090 kalmanFit->AddCalibration(transformation);
1091 //
1092 sprintf(tname,"tTPCTiltingZIROCOROCSideC_Sin%d",i);
1093 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1094 transformation->SetParams(0,0.03,0,&fpar);
1095 kalmanFit->AddCalibration(transformation);
1096 //
1097 //
1098 //
1099 fpar[2]=0; // absolute misalignment
1100 sprintf(tname,"tTPCTiltingZSectorSideA_Sin%d",i);
1101 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" , 0);
1102 transformation->SetParams(0,0.1,0,&fpar);
1103 kalmanFit->AddCalibration(transformation);
1104 //
1105 sprintf(tname,"tTPCTiltingZSectorSideC_Sin%d",i);
1106 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" , 0);
1107 transformation->SetParams(0,0.1,0,&fpar);
1108 kalmanFit->AddCalibration(transformation);
1109 }
171ba885 1110}
1111
1112
1113
1114void AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit){
1115 //
1116 //
1117 //
1118 TVectorD fpar(10);
1119 AliTPCTransformation * transformation=0;
1120 TBits maskInnerA(72);
1121 TBits maskInnerC(72);
1122 for (Int_t i=0;i<72;i++){
1123 if (i<36){
1124 if (i%36<18) maskInnerA[i]=kTRUE;
1125 if (i%36>=18) maskInnerC[i]=kTRUE;
1126 }
1127 }
1128 //
1129 //
1130 transformation = new AliTPCTransformation("tTPCDeltaLxIROCA", new TBits(maskInnerA), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1131 transformation->SetParams(0,0.2,0,&fpar);
1132 kalmanFit->AddCalibration(transformation);
1133 transformation = new AliTPCTransformation("tTPCDeltaLxIROCC", new TBits(maskInnerC), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1134 transformation->SetParams(0,0.2,0,&fpar);
1135 kalmanFit->AddCalibration(transformation);
1136 //
1137 transformation = new AliTPCTransformation("tTPCDeltaLyIROCA", new TBits(maskInnerA), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1138 transformation->SetParams(0,0.2,0,&fpar);
1139 kalmanFit->AddCalibration(transformation);
1140 transformation = new AliTPCTransformation("tTPCDeltaLyIROCC", new TBits(maskInnerC), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1141 transformation->SetParams(0,0.2,0,&fpar);
1142 kalmanFit->AddCalibration(transformation);
1143}
1144
d8b91e72 1145void AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit){
1146 //
1147 //
1148 //
1149 TVectorD fpar(10);
1150 AliTPCTransformation * transformation=0;
1151 Int_t fixSector =4;
1152 //
1153 for (Int_t isec=0; isec<36;isec++){
1154 TBits mask(72);
1155 mask[isec]=kTRUE;
1156 char tname[1000];
1157 //
1158 sprintf(tname,"tTPClocalLxIROC%d",isec);
1159 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1160 transformation->SetParams(0,0.2,0,&fpar);
1161 kalmanFit->AddCalibration(transformation);
1162 //
1163 sprintf(tname,"tTPClocalLyIROC%d",isec);
1164 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1165 transformation->SetParams(0,0.2,0,&fpar);
1166 kalmanFit->AddCalibration(transformation);
1167
1168 sprintf(tname,"tTPClocalRzIROC%d",isec);
1169 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0, 0);
1170 transformation->SetParams(0,0.2,0,&fpar);
1171 kalmanFit->AddCalibration(transformation);
1172
1173 }
1174 //
1175 for (Int_t isec=0; isec<36;isec++){
1176 if (isec%18==fixSector) continue;
1177 TBits mask(72);
1178 mask[isec] =kTRUE;
1179 mask[isec+36]=kTRUE;
1180 char tname[1000];
1181 //
1182 sprintf(tname,"tTPClocalLxSector%d",isec);
1183 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1184 transformation->SetParams(0,0.2,0,&fpar);
1185 kalmanFit->AddCalibration(transformation);
1186 //
1187 sprintf(tname,"tTPClocalLySector%d",isec);
1188 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1189 transformation->SetParams(0,0.2,0,&fpar);
1190 kalmanFit->AddCalibration(transformation);
1191
1192 sprintf(tname,"tTPClocalRzSector%d",isec);
1193 transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0, 0);
1194 transformation->SetParams(0,0.2,0,&fpar);
1195 kalmanFit->AddCalibration(transformation);
1196 }
1197
1198
1199}
1200
1201
1202void AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1203 //
1204 //
1205 //
1206 TVectorD fpar(10);
1207 AliTPCTransformation * transformation=0;
1208
1209 for (Int_t i=0; i<=ncos;i++){
1210 char tname[1000];
1211 fpar[0]=i; // cos frequency
1212 fpar[1]=0; // no sinus
1213 fpar[2]=1; // relative misalignment
1214 //
1215 // Local X shift
1216 //
1217 sprintf(tname,"tTPClocalLxIROCOROCSideA_Cos%d",i);
1218 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1219 transformation->SetParams(0,0.03,0,&fpar);
1220 kalmanFit->AddCalibration(transformation);
1221 sprintf(tname,"tTPClocalLxIROCOROCSideC_Cos%d",i);
1222 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1223 transformation->SetParams(0,0.03,0,&fpar);
1224 kalmanFit->AddCalibration(transformation);
1225 //
1226 //
1227 // Local Y shift
1228 //
1229 sprintf(tname,"tTPClocalLyIROCOROCSideA_Cos%d",i);
1230 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1231 transformation->SetParams(0,0.03,0,&fpar);
1232 kalmanFit->AddCalibration(transformation);
1233 sprintf(tname,"tTPClocalLyIROCOROCSideC_Cos%d",i);
1234 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1235 transformation->SetParams(0,0.03,0,&fpar);
1236 kalmanFit->AddCalibration(transformation);
1237 //
1238 //
1239 //
1240 // Z rotation
1241 //
1242 sprintf(tname,"tTPClocalRzIROCOROCSideA_Cos%d",i);
1243 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1244 transformation->SetParams(0,0.3,0,&fpar);
1245 kalmanFit->AddCalibration(transformation);
1246 sprintf(tname,"tTPClocalRzIROCOROCSideC_Cos%d",i);
1247 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1248 transformation->SetParams(0,0.3,0,&fpar);
1249 kalmanFit->AddCalibration(transformation);
1250 }
1251 //
1252 for (Int_t i=1; i<=nsin;i++){
1253 char tname[1000];
1254 fpar[0]=0; // cos frequency
1255 fpar[1]=i; // sinus frequency
1256 fpar[2]=1; // relative misalignment
1257 //
1258 // Local X shift
1259 //
1260 sprintf(tname,"tTPClocalLxIROCOROCSideA_Sin%d",i);
1261 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1262 transformation->SetParams(0,0.03,0,&fpar);
1263 kalmanFit->AddCalibration(transformation);
1264 sprintf(tname,"tTPClocalLxIROCOROCSideC_Sin%d",i);
1265 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1266 transformation->SetParams(0,0.03,0,&fpar);
1267 kalmanFit->AddCalibration(transformation);
1268 //
1269 //
1270 // Local Y shift
1271 //
1272 sprintf(tname,"tTPClocalLyIROCOROCSideA_Sin%d",i);
1273 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1274 transformation->SetParams(0,0.03,0,&fpar);
1275 kalmanFit->AddCalibration(transformation);
1276 sprintf(tname,"tTPClocalLyIROCOROCSideC_Sin%d",i);
1277 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1278 transformation->SetParams(0,0.03,0,&fpar);
1279 kalmanFit->AddCalibration(transformation);
1280 //
1281 //
1282 //
1283 // Z rotation
1284 //
1285 sprintf(tname,"tTPClocalRzIROCOROCSideA_Sin%d",i);
1286 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1287 transformation->SetParams(0,0.3,0,&fpar);
1288 kalmanFit->AddCalibration(transformation);
1289 sprintf(tname,"tTPClocalRzIROCOROCSideC_Sin%d",i);
1290 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1291 transformation->SetParams(0,0.3,0,&fpar);
1292 kalmanFit->AddCalibration(transformation);
1293 }
1294}
1295
1296void AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1297 //
1298 //
1299 //
1300 TVectorD fpar(10);
1301 AliTPCTransformation * transformation=0;
1302
1303 for (Int_t i=0; i<=ncos;i++){
1304 char tname[1000];
1305 fpar[0]=i; // cos frequency
1306 fpar[1]=0; // no sinus
1307 fpar[2]=0; // absolute misalignment
1308 if (i>0){
1309 //
1310 // A side is reference
1311 // local x
1312 sprintf(tname,"tTPClocalLxSectorSideA_Cos%d",i);
1313 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1314 transformation->SetParams(0,0.03,0,&fpar);
1315 kalmanFit->AddCalibration(transformation);
1316 if (i>1){
1317 //
1318 // Local Y shift
1319 //
1320 sprintf(tname,"tTPClocalLySectorSideA_Cos%d",i);
1321 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1322 transformation->SetParams(0,0.03,0,&fpar);
1323 kalmanFit->AddCalibration(transformation);
1324 }
1325 //
1326 //
1327 }
1328 //
1329 // C side to vary
1330 // local x
1331 sprintf(tname,"tTPClocalLxSectorSideC_Cos%d",i);
1332 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1333 transformation->SetParams(0,0.03,0,&fpar);
1334 kalmanFit->AddCalibration(transformation);
1335 //
1336 // Local Y shift
1337 //
1338 sprintf(tname,"tTPClocalLySectorSideC_Cos%d",i);
1339 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1340 transformation->SetParams(0,0.03,0,&fpar);
1341 kalmanFit->AddCalibration(transformation);
1342 //
1343 //
1344 // Z rotation - independent
1345 //
1346 sprintf(tname,"tTPClocalRzSectorSideA_Cos%d",i);
1347 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1348 transformation->SetParams(0,0.3,0,&fpar);
1349 kalmanFit->AddCalibration(transformation);
1350 sprintf(tname,"tTPClocalRzSectorSideC_Cos%d",i);
1351 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1352 transformation->SetParams(0,0.3,0,&fpar);
1353 kalmanFit->AddCalibration(transformation);
1354 }
1355
1356
1357
1358
1359 //
1360 //
1361 //
1362 for (Int_t i=1; i<=nsin;i++){
1363 char tname[1000];
1364 fpar[0]=0; // non cos frequency
1365 fpar[1]=i; // sinus frequency
1366 fpar[2]=0; // absolute misalignment
1367 //
1368 // Local X shift
1369 //
1370 sprintf(tname,"tTPClocalLxSectorSideA_Sin%d",i);
1371 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1372 transformation->SetParams(0,0.03,0,&fpar);
1373 kalmanFit->AddCalibration(transformation);
1374 sprintf(tname,"tTPClocalLxSectorSideC_Sin%d",i);
1375 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0, 0);
1376 transformation->SetParams(0,0.03,0,&fpar);
1377 kalmanFit->AddCalibration(transformation);
1378 //
1379 //
1380 // Local Y shift
1381 //
1382 sprintf(tname,"tTPClocalLySectorSideA_Sin%d",i);
1383 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1384 transformation->SetParams(0,0.03,0,&fpar);
1385 kalmanFit->AddCalibration(transformation);
1386 sprintf(tname,"tTPClocalLySectorSideC_Sin%d",i);
1387 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0, 0);
1388 transformation->SetParams(0,0.03,0,&fpar);
1389 kalmanFit->AddCalibration(transformation);
1390 //
1391 //
1392 //
1393 // Z rotation
1394 //
1395 sprintf(tname,"tTPClocalRzSectorSideA_Sin%d",i);
1396 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1397 transformation->SetParams(0,0.3,0,&fpar);
1398 kalmanFit->AddCalibration(transformation);
1399 sprintf(tname,"tTPClocalRzSectorSideC_Sin%d",i);
1400 transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0, 0);
1401 transformation->SetParams(0,0.3,0,&fpar);
1402 kalmanFit->AddCalibration(transformation);
1403 }
1404}
1405
1406void SelectPixel(){
1407 for (Int_t i=0;i<8;i++){
1408 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1409 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1410 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1411 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1412 }
1413}
1414
1415void SelectNonPixelA(){
1416 for (Int_t i=0;i<8;i++){
1417 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1418 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1419 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1420 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1421 }
1422}
1423
1424void SelectNonPixelC(){
1425 for (Int_t i=0;i<8;i++){
1426 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1427 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1428 kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1429 kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1430 }
1431}
1432
1433
1434void DumpQA1D( TObjArray &arrayOut){
1435 //
1436 //
1437 //
1438 TF1 fg("fg","gaus");
1439 TMatrixD sideARMS(8,2);
1440 TMatrixD sideCRMS(8,2);
1441 TMatrixD sideACRMS(8,2);
1442 TH1 *his=0;
1443 //
1444 // A side
1445 //
1446 SelectNonPixelA();
1447 for (Int_t i=0; i<8;i++){
1448 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1449 his->Fit(&fg,"","", -0.15,0.15);
1450 sideARMS(i,0) = fg.GetParameter(2);
1451 his->SetDirectory(0);
1452 his->SetName(Form("Original SideA_%s",his->GetName()));
1453 his->SetTitle(Form("Original SideA_%s",his->GetTitle()));
1454 arrayOut.AddLast(his);
1455 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1456 his->Fit(&fg,"","", -0.15,0.15);
1457 sideARMS(i,1) = fg.GetParameter(2);
1458 his->SetDirectory(0);
1459 his->SetName(Form("Aligned SideA_%s",his->GetName()));
1460 his->SetTitle(Form("Aligned SideA_%s",his->GetTitle()));
1461 arrayOut.AddLast(his);
1462 }
1463 //
1464 // C side
1465 //
1466 SelectNonPixelC();
1467 for (Int_t i=0; i<8;i++){
1468 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1469 his->Fit(&fg,"","", -0.15,0.15);
1470 sideCRMS(i,0) = fg.GetParameter(2);
1471 his->SetDirectory(0);
1472 his->SetName(Form("Original SideC_%s",his->GetName()));
1473 his->SetTitle(Form("Original SideC_%s",his->GetTitle()));
1474 arrayOut.AddLast(his);
1475 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1476 his->Fit(&fg,"","", -0.15,0.15);
1477 sideCRMS(i,1) = fg.GetParameter(2);
1478 his->SetDirectory(0);
1479 his->SetName(Form("Aligned SideC_%s",his->GetName()));
1480 his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1481 arrayOut.AddLast(his);
1482 }
1483 //
1484 // AC side
1485 //
1486 SelectPixel();
1487 for (Int_t i=0; i<8;i++){
1488 his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1489 his->Fit(&fg,"","", -0.15,0.15);
1490 sideACRMS(i,0) = fg.GetParameter(2);
1491 his->SetDirectory(0);
1492 his->SetName(Form("Original SideAC_%s",his->GetName()));
1493 his->SetTitle(Form("Original SideAC_%s",his->GetTitle()));
1494 arrayOut.AddLast(his);
1495 his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1496 his->Fit(&fg,"","", -0.15,0.15);
1497 sideACRMS(i,1) = fg.GetParameter(2);
1498 his->SetDirectory(0);
1499 his->SetName(Form("Aligned SideC_%s",his->GetName()));
1500 his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1501 arrayOut.AddLast(his);
1502 }
1503 printf("DumQA\n");
1504 sideARMS.Print();
1505 sideCRMS.Print();
1506 sideACRMS.Print();
1507}
1508
1509void MakeFolders(TObjArray * arrayOut){
1510 //
1511 //
1512 //
1513 TFolder *folderBase = new TFolder("TPC align","TPC align");
1514 //
1515 //
1516 TString maskDelta[8];
1517 TString maskAlign[2]={"Orig","Alig");
1518 for (Int_t i=0; i<8;i++){
1519 maskDelta[i] = kalmanFitNew->fLinearTrackDelta[i]->GetName();
1520 }
1521
1522 Int_t entries = arrayOut->GetEntries();
1523 //
1524}
1525
1526
1527
1528
171ba885 1529void MergeKalman(const char * list = "kalmanFit.list"){
1530 //
1531 //
1532 //
1533 ifstream in;
1534 in.open(list);
1535 TString currentFile;
d8b91e72 1536 kalmanFitNew= 0;
171ba885 1537 Int_t counter=0;
1538 while(in.good()) {
d8b91e72 1539 //
1540 // calibrated
1541 //
1542 in >> currentFile;
171ba885 1543 printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1544 if (currentFile.Length()==0) continue;
1545 TFile * ffit = TFile::Open(currentFile.Data());
d8b91e72 1546 TEntryList * tlist = (TEntryList*) ffit->Get("eventList");
1547 TMatrixD * cMatrix = (TMatrixD*) ffit->Get("cutMarix");
1548 if (tlist&&cMatrix){
1549 printf("Track entries=%d\n",tlist->GetN());
1550 if (cMatrix->Sum()<=0.00000001) {
1551 printf("Problem with track selection\n");
1552 continue;
1553 }
1554 }else{
1555 printf("Problem with track selection\n");
1556 continue;
1557 }
1558 //
1559
171ba885 1560 AliTPCkalmanFit * fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFit");
1561 if (!fit) continue;
d8b91e72 1562 if (!kalmanFitNew) {kalmanFitNew= fit; continue;};
1563 kalmanFitNew->Add(fit);
1564 printf("Selected entries=%f\n",fit->fLinearTrackDelta[0]->GetEntries());
1565 //delete tlist;
1566 //delete cMatrix;
1567 delete fit;
1568 delete ffit;
1569 //
1570 // original
1571 //
1572 currentFile.ReplaceAll("TPC","TPCOrig");
1573 printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1574 if (currentFile.Length()==0) continue;
1575 ffit = TFile::Open(currentFile.Data());
1576 fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFitOrig");
1577 if (!fit) continue;
1578 if (!kalmanFitOrig) {kalmanFitOrig= fit; continue;};
1579 kalmanFitOrig->Add(fit);
171ba885 1580 delete fit;
1581 delete ffit;
1582 counter++;
1583 }
d8b91e72 1584 //
1585 // save merged results
1586 //
1587 TFile f("mergeKalmanFit.root","recreate");
1588 kalmanFitNew->Write("kalmanFitNew");
1589 kalmanFitOrig->Write("kalmanFitOrig");
1590 f.Close();
171ba885 1591}
1592
1593
d8b91e72 1594
171ba885 1595/*
d8b91e72 1596 Example - submit alignment as batch jobs
1597 ifile=0;
1598 ntracksSkip=200
1599 ntracksSkipOffset=0
1600 nclustersSkip=2
1601 nclustersSkipOffset=0
1602 ntracks=1000000
1603 ndump=5
1604 isTest=0
171ba885 1605 bDir=`pwd`
d8b91e72 1606 ver=aaa;
1607
1608 while [ $ntracksSkipOffset -lt $ntracksSkip ] ; do
1609 nclustersSkipOffset=0;
1610 while [ $nclustersSkipOffset -lt $nclustersSkip ] ; do
1611 echo Tr $ntracksSkipOffset Cl $nclustersSkipOffset;
1612 ver=kalmanDirTrack$ntracksSkipOffset$nclustersSkipOffset
1613 echo New Directory $ver
1614 mkdir $ver;
1615 cd $ver;
1616 cp $bDir/align.txt . ;
1617 ln -sf $bDir/kalmanFitApply.root . ;
1618 echo aliroot -q -b "$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C($ntracks,10000,0,$ndump,$ntracksSkip,$ntracksSkipOffset, $nclustersSkip,$nclustersSkipOffset,$isTest)" ;
1619 bsub -q alice-t3_8h -o `pwd`/output.log command aliroot -q -b "$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C($ntracks,10000,0,$ndump,$ntracksSkip,$ntracksSkipOffset, $nclustersSkip,$nclustersSkipOffset,$isTest)" ;
1620 nclustersSkipOffset=$(( $nclustersSkipOffset + 1 ))
1621 cd $bDir;
1622 echo $bDir;
1623 done;
1624 cd $bDir;
1625 echo $bDir;
1626 ntracksSkipOffset=$(( $ntracksSkipOffset + 1 ))
1627 echo Tr $ntracksSkipOffset;
1628done
171ba885 1629
171ba885 1630
1631
d8b91e72 1632*/
1633
171ba885 1634
1635
1636
1637
1638
1639
1640
1641
1642
1643