]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCkalmanAlign.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bastian's Task for Jet Fragmentation, to be merged with other FF tasks
[u/mrichter/AliRoot.git] / TPC / AliTPCkalmanAlign.cxx
CommitLineData
125d3a38 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17 TPC Kalman filter implementation for TPC sector alignment.
18 Output of the AliTPCcalibAlign is used as a input for TPC global alignment.
19 In AliTPCcalibAlign histograms - track parameter matching at sector boundaries are created.
20 Each sector is aligned with 5 neighborhoud (sectors)
21 1. Up-down - 1
22 2. Left-right - 4
23
24 Sector alignment parameters are obtained finding the alignment parameters, minimizing
25 misalignmnet for all piars fo sectors.
26
27 Global minimization- MakeGlobalAlign
28
29
30 Example usage:
31 gSystem->Load("libANALYSIS");
32 gSystem->Load("libTPCcalib");
33 //
4486a91f 34 Int_t run=117220;
125d3a38 35 .x ConfigCalibTrain.C(run)
4486a91f 36 calibDB = AliTPCcalibDB::Instance()
125d3a38 37
38 AliTPCkalmanAlign kalmanAlign("TPC align", "TPC align"); // create the object
39 kalmanAlign.ReadAlign("CalibObjects.root"); // read the calibration form file
40 kalmanAlign.MakeGlobalAlign(); // do kalman alignment
41 kalmanAlign.DrawDeltaAlign(); // make QA plot
42 //
43
44
45*/
46#include "TMath.h"
47#include "TTreeStream.h"
48#include "TGraph.h"
49#include "TGraphErrors.h"
50#include "TVectorD.h"
51#include "TClonesArray.h"
4486a91f 52#include "TCut.h"
53#include "TChain.h"
54#include "AliXRDPROOFtoolkit.h"
55#include "TLegend.h"
56#include "TCanvas.h"
125d3a38 57
82628455 58#include "AliTPCcalibDB.h"
59#include "AliTPCCalROC.h"
125d3a38 60#include "AliCDBMetaData.h"
61#include "AliCDBId.h"
62#include "AliCDBManager.h"
63#include "AliCDBStorage.h"
64#include "AliCDBEntry.h"
65#include "AliAlignObjParams.h"
66#include "AliTPCROC.h"
82628455 67#include "AliTracker.h"
125d3a38 68#include "TFile.h"
69#include "TLinearFitter.h"
70#include "AliTPCcalibAlign.h"
71#include "TH1.h"
72#include "AliTPCCalPad.h"
73#include "AliTPCkalmanAlign.h"
4486a91f 74#include "TStatToolkit.h"
75#include "AliTPCPreprocessorOnline.h"
76#include "TPostScript.h"
125d3a38 77
78AliTPCkalmanAlign::AliTPCkalmanAlign():
79 TNamed(),
80 fCalibAlign(0), // kalman alignnmnt
81 fOriginalAlign(0), // original alignment 0 read for the OCDB
4486a91f 82 fNewAlign(0),
83 fPadTime0(0),
84 fFitCEGlobal(0),
85 fFitCELocal(0)
125d3a38 86{
87 //
88 // Default constructor
89 //
90 for (Int_t i=0; i<4; i++){
91 fDelta1D[i]=0;
92 fCovar1D[i]=0;
93 }
94}
95
96AliTPCkalmanAlign::AliTPCkalmanAlign(const char* name, const char* title):
97 TNamed(name,title),
98 fCalibAlign(0), // kalman alignnmnt
99 fOriginalAlign(0), // original alignment 0 read for the OCDB
4486a91f 100 fNewAlign(0),
101 fPadTime0(0),
102 fFitCEGlobal(0),
103 fFitCELocal(0)
125d3a38 104{
105 //
106 // Default constructor
107 //
108 for (Int_t i=0; i<4; i++){
109 fDelta1D[i]=0;
110 fCovar1D[i]=0;
111 }
4486a91f 112 fFitCEGlobal = new TObjArray(6);
113 fFitCELocal = new TObjArray(6);
114 for (Int_t ipar=0; ipar<6;ipar++){
115 fFitCEGlobal->AddAt(new TVectorD(36),ipar);
116 fFitCELocal->AddAt(new TVectorD(36),ipar);
117 }
125d3a38 118}
119
120void AliTPCkalmanAlign::ReadAlign(const char *fname){
121 //
122 // Read old alignment used in the reco
123 // and the residual histograms
124 // WE ASSUME that the OCDB path is set in the same way as done in the calibration
125 //
126 TFile fcalib(fname);
127 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
128 fCalibAlign=0;
129 if (array) fCalibAlign=( AliTPCcalibAlign *)array->FindObject("alignTPC");
130 fCalibAlign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
131 //
132 // old alignment used
133 AliCDBEntry * cdbEntry= AliCDBManager::Instance()->Get("TPC/Align/Data");
134 fOriginalAlign =0;
135 if (cdbEntry){
136 fOriginalAlign = (TClonesArray*)cdbEntry->GetObject();
137 }
138
139}
140
141void AliTPCkalmanAlign::BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t mean, Double_t sigma){
142 //
143 // Book Align 1D parameters and covariance
144 //
145 param.ResizeTo(72,1);
146 covar.ResizeTo(72,72);
147 for (Int_t i=0;i<72;i++){
148 param(i,0)=mean;
149 for (Int_t j=0;j<72;j++) covar(i,j)=0;
150 covar(i,i)=sigma*sigma;
151 }
152}
153
154
155void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &vecXk, TMatrixD &covXk){
156 //
157 // Update 1D kalman filter
158 //
159 const Int_t knMeas=1;
160 const Int_t knElem=72;
161 TMatrixD mat1(72,72); // update covariance matrix
162 TMatrixD matHk(1,knElem); // vector to mesurement
163 TMatrixD vecYk(knMeas,1); // Innovation or measurement residual
164 TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose
165 TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance
166 TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain
167 TMatrixD covXk2(knElem,knElem); // helper matrix
168 TMatrixD covXk3(knElem,knElem); // helper matrix
169 TMatrixD vecZk(1,1);
170 TMatrixD measR(1,1);
171 vecZk(0,0)=delta;
172 measR(0,0)=sigma*sigma;
173 //
174 // reset matHk
175 for (Int_t iel=0;iel<knElem;iel++)
176 for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
177 //mat1
178 for (Int_t iel=0;iel<knElem;iel++) {
179 for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
180 mat1(iel,iel)=1;
181 }
182 //
183 matHk(0, s1)=1;
184 matHk(0, s2)=-1;
185 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
186 matHkT=matHk.T(); matHk.T();
187 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
188 matSk.Invert();
189 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
190 vecXk += matKk*vecYk; // updated vector
191 covXk2= (mat1-(matKk*matHk));
192 covXk3 = covXk2*covXk;
193 covXk = covXk3;
194}
195
196
197
198
199void AliTPCkalmanAlign::MakeGlobalAlign(){
200 //
201 // Combine all pairs of fitters and make global alignemnt
202 //
203 AliTPCkalmanAlign &kalmanAlign=*this;
204 TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
82628455 205 //
206 // get ce info
207 //
208 AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
209 TVectorD paramCE[72];
210 TMatrixD covarCE[72];
4486a91f 211 Int_t statUpDown=0; // statistic up down
212 Int_t statLeftRight=0; // statistic left-right
82628455 213 Float_t chi2;
214 for (Int_t isec=0; isec<72; isec++){
215 AliTPCCalROC * roc0 = padTime0->GetCalROC(isec);
216 roc0->GlobalFit(0,kFALSE,paramCE[isec],covarCE[isec],chi2,0);
217 (*pcstream)<<"ceFit"<<
218 "isec="<<isec<<
219 "p0.="<<&paramCE[isec]<<
220 "\n";
221 }
222
125d3a38 223 DumpOldAlignment(pcstream);
224 const Int_t kMinEntries=400;
225 TMatrixD vec[5];
226 TMatrixD cov[5];
227 kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
228 kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
229 kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
230 kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
231 //
232 TVectorD delta(5);
233 TVectorD rms(5);
82628455 234 TVectorD stat(5);
125d3a38 235 TH1 * his=0;
236 for (Int_t is0=0;is0<72;is0++)
237 for (Int_t is1=0;is1<72;is1++){
238 //
239 //
240 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
241 if (!his) continue;
242 if (his->GetEntries()<kMinEntries) continue;
243 delta[0]=his->GetMean();
82628455 244 rms[0]=his->GetRMS();
245 stat[0]=his->GetEntries();
125d3a38 246 kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
247 //
248 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
249 if (!his) continue;
250 delta[1]=his->GetMean();
82628455 251 rms[1]=his->GetRMS();
252 stat[1]=his->GetEntries();
125d3a38 253 kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
254 //
255 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
256 if (!his) continue;
257 delta[2] = his->GetMean();
82628455 258 rms[2]=his->GetRMS();
259 stat[2]=his->GetEntries();
125d3a38 260 kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
261 //
262 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
263 if (!his) continue;
264 delta[3] = his->GetMean();
82628455 265 rms[3]=his->GetRMS();
266 stat[3]=his->GetEntries();
125d3a38 267 kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
4486a91f 268 if (is1==is0+36) statUpDown+=Int_t(stat[0]);
269 if (is1==is0+35) statLeftRight+=Int_t(stat[0]);
125d3a38 270 }
82628455 271
272 fDelta1D[0] = (TMatrixD*)vec[0].Clone();
273 fDelta1D[1] = (TMatrixD*)vec[1].Clone();
274 fDelta1D[2] = (TMatrixD*)vec[2].Clone();
275 fDelta1D[3] = (TMatrixD*)vec[3].Clone();
276 //
277 fCovar1D[0] = (TMatrixD*)cov[0].Clone();
278 fCovar1D[1] = (TMatrixD*)cov[1].Clone();
279 fCovar1D[2] = (TMatrixD*)cov[2].Clone();
280 fCovar1D[3] = (TMatrixD*)cov[3].Clone();
4486a91f 281 statUpDown/=36;
282 statLeftRight/=36;
82628455 283 MakeNewAlignment(kTRUE);
4486a91f 284 //FitCE();
125d3a38 285 for (Int_t is0=0;is0<72;is0++)
286 for (Int_t is1=0;is1<72;is1++){
82628455 287 Bool_t isPair=kFALSE;
288 if (TMath::Abs(is0%18-is1%18)<2) isPair=kTRUE;
289 if (TMath::Abs(is0%18-is1%18)==17) isPair=kTRUE;
290 if (!isPair) continue;
291 stat[0]=0; stat[1]=0; stat[2]=0; stat[3]=0;
125d3a38 292 //
293 //
294 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
82628455 295 if (his){
296 delta[0]=his->GetMean();
297 rms[0]=his->GetRMS();
298 stat[0]=his->GetEntries();
299 }
125d3a38 300 //
301 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
82628455 302 if (his) {
303 delta[1]=his->GetMean();
304 rms[1]=his->GetRMS();
305 stat[1]=his->GetEntries();
306 }
125d3a38 307 //
308 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
82628455 309 if (his){
310 delta[2] = his->GetMean();
311 rms[2]=his->GetRMS();
312 stat[2]=his->GetEntries();
313 }
125d3a38 314 //
315 his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
82628455 316 if (his){
317 delta[3] = his->GetMean();
318 rms[3]=his->GetRMS();
319 stat[3]=his->GetEntries();
320 }
321 Int_t run = AliCDBManager::Instance()->GetRun();
322 Float_t bz = AliTracker::GetBz();
4486a91f 323 TVectorD fceG[6],fceL[6];
324 for (Int_t ipar=0; ipar<6;ipar++){
325 fceG[ipar]=*((TVectorD*)fFitCEGlobal->At(ipar));
326 fceL[ipar]=*((TVectorD*)fFitCELocal->At(ipar));
327 }
125d3a38 328 (*pcstream)<<"kalmanAlignDebug"<<
82628455 329 "run="<<run<<
330 "bz="<<bz<<
125d3a38 331 "is0="<<is0<<
332 "is1="<<is1<<
333 "delta.="<<&delta<<
82628455 334 "rms.="<<&rms<<
335 "stat.="<<&stat<<
125d3a38 336 "vec0.="<<&vec[0]<<
337 "vec1.="<<&vec[1]<<
338 "vec2.="<<&vec[2]<<
339 "vec3.="<<&vec[3]<<
82628455 340 "pceIn0.="<<&paramCE[is0%36]<<
341 "pceOut0.="<<&paramCE[is0%36+36]<<
342 "pceIn1.="<<&paramCE[is1%36]<<
343 "pceOut1.="<<&paramCE[is1%36+36]<<
4486a91f 344 "fceG0.="<<&fceG[0]<< // global fit of CE
345 "fceG1.="<<&fceG[1]<< // global fit of CE
346 "fceG2.="<<&fceG[2]<< // global fit of CE
347 "fceG3.="<<&fceG[3]<< // global fit of CE
348 "fceG4.="<<&fceG[4]<< // global fit of CE
349 "fceL5.="<<&fceG[5]<< // global fit of CE
350 "fceL0.="<<&fceL[0]<< // global fit of CE
351 "fceL1.="<<&fceL[1]<< // global fit of CE
352 "fceL2.="<<&fceL[2]<< // global fit of CE
353 "fceL3.="<<&fceL[3]<< // global fit of CE
354 "fceL4.="<<&fceL[4]<< // global fit of CE
355 "fceL5.="<<&fceL[5]<< // global fit of CE
125d3a38 356 "\n";
357 }
4486a91f 358
359 Int_t run = AliCDBManager::Instance()->GetRun();
360 Float_t bz = AliTracker::GetBz();
361 (*pcstream)<<"runSummary"<<
362 "run="<<run<< // run number
363 "bz="<<bz<< // bz field
364 "statUpDown="<<statUpDown<< // stat up-down
365 "statLeftRight="<<statLeftRight<< // stat left-right
366 "\n";
367
125d3a38 368 delete pcstream;
369}
370
371
372
373
374
375
376void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad *pad, Int_t ustartRun, Int_t uendRun, const char* storagePath ){
377 //
378 // Update OCDB
379 // .x ConfigCalibTrain.C(117116)
380 // AliTPCcalibDB::Instance()->GetPulserTmean()
381 // pad->Add(-pad->GetMean())
382 AliCDBMetaData *metaData= new AliCDBMetaData();
383 metaData->SetObjectClassName("TObjArray");
384 metaData->SetResponsible("Marian Ivanov");
385 metaData->SetBeamPeriod(1);
386 metaData->SetAliRootVersion("05-25-01"); //root version
387 metaData->SetComment("Calibration of the z - time misalignment");
388 AliCDBId* id1=NULL;
389 id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
390 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
391 gStorage->Put(pad, (*id1), metaData);
392}
393
394
395
396void AliTPCkalmanAlign::DrawDeltaAlign(){
397 //
398 // Draw the reuslts of the alignment
399 // Residual misalignment in respect with previous alignment shown
400 //
401 //
402 TFile f("AliTPCkalmanAlign.root","update");
403 TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
404 TH1::AddDirectory(0);
405 //
406 treeDelta->SetAlias("sector","is0");
407 treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
408 treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
409 treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
410 treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
411 //
412 treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
413 treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
414 treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
415 treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
416 //
417 treeDelta->SetMarkerStyle(25);
418 treeDelta->SetMarkerColor(4);
419 treeDelta->SetLineColor(4);
420 const char *type[3]={"up-down","left-right","right-left"};
421 const char *gropt[3]={"alp","lp","lp"};
422 const char *varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10*vec0.fElements[is0]:sector"};
423 const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
424 const char *varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10*vec1.fElements[is0]:sector"};
425 const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
426 const char *varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000*vec3.fElements[is0]:sector"};
427 const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
428 const char *varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000*vec2.fElements[is0]:sector"};
429 const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
430 TLegend *legend = 0;
431 Int_t grcol[3]={2,1,4};
432 Int_t entries=0;
433 TGraph *grDelta[3]={0,0,0};
434 TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
435 TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
436 TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
437 TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
438 canvasDy->Divide(2,2);
439 canvasDz->Divide(2,2);
440 canvasDtheta->Divide(2,2);
441 canvasDphi->Divide(2,2);
442
443 //
444 // Dy
445 //
446 canvasDy->cd(1);
447 treeDelta->Draw("dYmeas:dYfit");
448 for (Int_t itype=0; itype<3; itype++){
449 canvasDy->cd(itype+2);
450 entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
451 grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
452 entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
453 grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
454 entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
455 grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
456 legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
457 for (Int_t i=0; i<3; i++) {
458 grDelta[i]->SetMaximum(1.5);
459 grDelta[i]->SetMinimum(-1);
460 grDelta[i]->SetTitle(type[i]);
461 grDelta[i]->SetMarkerColor(grcol[i]);
462 grDelta[i]->SetLineColor(grcol[i]);
463 grDelta[i]->SetMarkerStyle(25+i);
464 grDelta[i]->GetXaxis()->SetTitle("sector");
465 grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)");
466 if (itype==2 && i>0) continue;
467 grDelta[i]->Draw(gropt[i]);
468 legend->AddEntry(grDelta[i]);
469 }
470 legend->Draw();
471 }
472 //
473 // Dz
474 //
475 canvasDz->cd();
476 canvasDz->cd(1);
477 treeDelta->Draw("dZmeas:dZfit");
478 for (Int_t itype=0; itype<3; itype++){
479 canvasDz->cd(itype+2);
480 entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
481 grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
482 entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
483 grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
484 entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
485 grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
486 legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
487 for (Int_t i=0; i<3; i++) {
488 grDelta[i]->SetMaximum(1.5);
489 grDelta[i]->SetMinimum(-1);
490 grDelta[i]->SetTitle(type[i]);
491 grDelta[i]->SetMarkerColor(grcol[i]);
492 grDelta[i]->SetLineColor(grcol[i]);
493 grDelta[i]->SetMarkerStyle(25+i);
494 grDelta[i]->GetXaxis()->SetTitle("sector");
495 grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)");
496 if (itype==2 && i>0) continue;
497 grDelta[i]->Draw(gropt[i]);
498 legend->AddEntry(grDelta[i]);
499 }
500 legend->Draw();
501 }
502
503 //
504 // Dtheta
505 //
506 canvasDtheta->cd(1);
507 treeDelta->Draw("dThetameas:dThetafit");
508 for (Int_t itype=0; itype<3; itype++){
509 canvasDtheta->cd(itype+2);
510 entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
511 grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
512 entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
513 grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
514 entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
515 grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
516 legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
517 for (Int_t i=0; i<3; i++) {
518 grDelta[i]->SetMaximum(4.);
519 grDelta[i]->SetMinimum(-3.);
520 grDelta[i]->SetTitle(type[i]);
521 grDelta[i]->SetMarkerColor(grcol[i]);
522 grDelta[i]->SetLineColor(grcol[i]);
523 grDelta[i]->SetMarkerStyle(25+i);
524 grDelta[i]->GetXaxis()->SetTitle("sector");
525 grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)");
526 if (itype==2 && i>0) continue;
527 grDelta[i]->Draw(gropt[i]);
528 legend->AddEntry(grDelta[i]);
529 }
530 legend->Draw();
531 }
532
533 //
534 // Dphi
535 //
536 canvasDphi->cd();
537 canvasDphi->cd(1);
538 treeDelta->Draw("dPhimeas:dPhifit");
539 for (Int_t itype=0; itype<3; itype++){
540 canvasDphi->cd(itype+2);
541 entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
542 grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
543 entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
544 grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
545 entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
546 grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
547 legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
548 for (Int_t i=0; i<3; i++) {
549 grDelta[i]->SetMaximum(4.);
550 grDelta[i]->SetMinimum(-3.);
551 grDelta[i]->SetTitle(type[i]);
552 grDelta[i]->SetMarkerColor(grcol[i]);
553 grDelta[i]->SetLineColor(grcol[i]);
554 grDelta[i]->SetMarkerStyle(25+i);
555 grDelta[i]->GetXaxis()->SetTitle("sector");
556 grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
557 if (itype==2 && i>0) continue;
558 grDelta[i]->Draw(gropt[i]);
559 legend->AddEntry(grDelta[i]);
560 }
561 legend->Draw();
562 }
563 //
564 //
565 f.cd();
566 canvasDy->Write();
567 canvasDz->Write();
568 canvasDtheta->Write();
569 canvasDphi->Write();
4486a91f 570 //
571 //
572 //
573 TPostScript *ps = new TPostScript("alignReport.ps", 112);
574 ps->NewPage();
575 canvasDy->Draw();
576 ps->NewPage();
577 canvasDz->Draw();
578 ps->NewPage();
579 canvasDtheta->Draw();
580 ps->NewPage();
581 canvasDphi->Draw();
582 ps->Close();
583 delete ps;
125d3a38 584}
585
586
587
588void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
589 // Dump the content of old alignemnt
590 // Expected that the old alignmnet is loaded
591 //
592 if (!fOriginalAlign) return;
593 //
594 TVectorD localTrans(3);
595 TVectorD globalTrans(3);
596 TVectorD localRot(3);
597 TVectorD globalRot(3);
598 AliGeomManager::ELayerID idLayer;
599 Int_t idModule=0;
600 //
601 for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
602 AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
603 params->GetVolUID(idLayer,idModule);
604 params->GetLocalTranslation(localTrans.GetMatrixArray());
605 params->GetLocalAngles(localRot.GetMatrixArray());
606 params->GetTranslation(globalTrans.GetMatrixArray());
607 params->GetAngles(globalRot.GetMatrixArray());
608 Int_t sector=idModule;
609 if (idLayer>7) sector+=36;
610 (*pcstream)<<"oldAlign"<<
611 //"idLayer="<<idLayer<<
612 "idModule="<<idModule<<
613 "sector="<<sector<<
614 "lT.="<<&localTrans<<
615 "gT.="<<&localTrans<<
616 "lR.="<<&localRot<<
617 "gR.="<<&globalRot<<
618 "\n";
619 }
620}
621
622
623void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
624 //
625 // make a new Alignment entry
626 //
627 if (!fOriginalAlign) return;
628 //
629 TVectorD localTrans(3);
630 TVectorD globalTrans(3);
631 TVectorD localRot(3);
632 TVectorD globalRot(3);
633 //
634 TVectorD localTransNew(3); // new entries
635 TVectorD globalTransNew(3);
636 TVectorD localRotNew(3);
637 TVectorD globalRotNew(3);
638 //
639 AliGeomManager::ELayerID idLayer;
640 Int_t idModule=0;
641 //
642 fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
643 for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
644 AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
82628455 645 //AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
125d3a38 646 params->GetVolUID(idLayer,idModule);
647 Int_t sector=(Int_t)idModule;
648 if (idLayer>7) sector+=36;
649 params->GetLocalTranslation(localTrans.GetMatrixArray());
650 params->GetLocalAngles(localRot.GetMatrixArray());
651 params->GetTranslation(globalTrans.GetMatrixArray());
652 params->GetAngles(globalRot.GetMatrixArray());
653 //
654 //
655 //
656 if (badd){ // addition if
657 localTransNew=localTrans;
658 localRotNew=localRot;
659 }
82628455 660 localTransNew[1]=localTransNew[1]-((*fDelta1D[0])(sector,0));
661 localRot[0] =localRot[0]-(*fDelta1D[2])(sector,0);
125d3a38 662 //
82628455 663 if (pcstream) (*pcstream)<<"alignParams"<<
125d3a38 664 //"idLayer="<<idLayer<<
665 "idModule="<<idModule<<
666 "sector="<<sector<<
667 "olT.="<<&localTrans<<
125d3a38 668 "olR.="<<&localRot<<
82628455 669 "ogT.="<<&localTrans<<
125d3a38 670 "ogR.="<<&globalRot<<
82628455 671 "nlT.="<<&localTransNew<<
672 "nlR.="<<&localRotNew<<
673 "ngT.="<<&localTransNew<<
674 "ngR.="<<&globalRotNew<<
125d3a38 675 "\n";
676 }
4486a91f 677}
678
679
680
681void AliTPCkalmanAlign::DrawAlignmentTrends(){
682 //
683 // Draw trends of alingment variables
684 //
685 /*
686 */
687 AliXRDPROOFtoolkit toolkit;
688 TChain * chain = toolkit.MakeChainRandom("align.list.Good","kalmanAlignDebug",0,2000);
689 TChain * chainRef = toolkit.MakeChainRandom("alignRef.list","kalmanAlignDebug",0,2000);
690 chain->AddFriend(chainRef,"R");
691 chainRef->AddFriend(chainRef,"T");
692 //cuts
693 TCut cutS="stat.fElements[0]>200&&stat.fElements[1]>200&&stat.fElements[3]>200&&stat.fElements[3]>200"; //statistic in the bin
694 TCut cutST="T.stat.fElements[0]>200&&T.stat.fElements[1]>200&&T.stat.fElements[3]>200&&T.stat.fElements[3]>200"; //statistic in the bin
695 // TTree *tree = chain->CopyTree(cutS);
696 //TTree *treeR = chainRef->CopyTree(cutST);
697
698 TCanvas * canvasDy= new TCanvas("canvasDy","canvasDy");
699 TH1 *his=0;
700 TLegend *legend = 0;
701 // Int_t grcol[3]={2,1,4};
125d3a38 702
4486a91f 703 legend = new TLegend(0.7,0.6,0.9,0.9, "Alignment #Delta_{y}- Up-Down");
704 for (Int_t isec=0; isec<18; isec+=2){
705 chain->SetMarkerColor(1+(isec%5));
706 chain->SetMarkerStyle(isec+20);
707 chain->Draw("10*(delta.fElements[0]-R.delta.fElements[0]):run",cutS+Form("is1==is0+36&&is0==%d",isec),"profgoff");
708 his = (TH1*)(chain->GetHistogram()->Clone());
709 his->SetName(Form("#Delta_{Y} sector %d",isec));
710 his->SetTitle(Form("#Delta_{Y} sector %d",isec));
711 his->SetMaximum(1.);
712 his->SetMinimum(-1.);
713 his->GetYaxis()->SetTitle("#Delta_{y} (mm)");
714 his->GetXaxis()->SetTitle("run Number");
715 if (isec==0) his->Draw("");
716 if (isec>0) his->Draw("same");
717 legend->AddEntry(his);
718 }
719 legend->Draw();
720 canvasDy->Draw();
721}
722
723
724
725
726
727
728void AliTPCkalmanAlign::FitCE(){
729 //
730 // fit CE
731 // 1. Global fit - gy and gx
732 // 2. Local X fit common
733 // 3. Sector fit
734 //
735 AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
736 //
737 AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
738 AliTPCCalPad *padNoise = AliTPCcalibDB::Instance()->GetPadNoise();
739 AliTPCCalPad * ceTmean = AliTPCcalibDB::Instance()->GetCETmean(); // CE information
740 AliTPCCalPad * ceTrms = AliTPCcalibDB::Instance()->GetCETrms();
741 AliTPCCalPad * ceQmean = AliTPCcalibDB::Instance()->GetCEQmean();
742 AliTPCCalPad * pulserTmean = AliTPCcalibDB::Instance()->GetPulserTmean(); //
743 AliTPCCalPad * pulserTrms = AliTPCcalibDB::Instance()->GetPulserTrms();
744 AliTPCCalPad * pulserQmean = AliTPCcalibDB::Instance()->GetPulserQmean();
745 AliTPCCalPad * dmap0 = AliTPCcalibDB::Instance()->GetDistortionMap(0); // distortion maps
746 AliTPCCalPad * dmap1 = AliTPCcalibDB::Instance()->GetDistortionMap(1);
747 AliTPCCalPad * dmap2 = AliTPCcalibDB::Instance()->GetDistortionMap(2);
748 pulserTmean->Add(-pulserTmean->GetMean());
749 //
750 preprocesor->AddComponent(padTime0->Clone());
751 preprocesor->AddComponent(padNoise->Clone());
752 preprocesor->AddComponent(pulserTmean->Clone());
753 preprocesor->AddComponent(pulserQmean->Clone());
754 preprocesor->AddComponent(pulserTrms->Clone());
755 preprocesor->AddComponent(ceTmean->Clone());
756 preprocesor->AddComponent(ceQmean->Clone());
757 preprocesor->AddComponent(ceTrms->Clone());
758 preprocesor->AddComponent(dmap0->Clone());
759 preprocesor->AddComponent(dmap1->Clone());
760 preprocesor->AddComponent(dmap2->Clone());
761 preprocesor->DumpToFile("cetmean.root");
762
763 TCut cutNoise="abs(PadNoise.fElements/PadNoise_Median-1)<0.3";
764 TCut cutPulserT="abs(PulserTrms.fElements/PulserTrms_Median-1)<0.2";
765 TCut cutPulserQ="abs(PulserQmean.fElements/PulserQmean_Median-1)<0.2";
766 TCut cutCEQ="CEQmean.fElements>50";
767 TCut cutCET="abs(CETmean.fElements)<2";
768 TCut cutAll=cutNoise+cutPulserT+cutPulserQ+cutCEQ+cutCET;
769 //
770 //
771 TFile * f = new TFile("cetmean.root");
772 TTree * chain = (TTree*) f->Get("calPads");
773 Int_t entries = chain->Draw("1",cutAll,"goff");
774 if (entries<200000) return; // no calibration available - pulser or CE or noise
775
776 TStatToolkit toolkit;
777 Double_t chi2=0;
778 Int_t npoints=0;
779 TVectorD param;
780 TMatrixD covar;
781 //
782 // make a aliases
783 AliTPCkalmanAlign::MakeAliasCE(chain);
784 TString fstringG=""; // global part
785 //
786 fstringG+="Gy++"; // par 1 - global y
787 fstringG+="Gx++"; // par 2 - global x
788 //
789 fstringG+="isin++"; // delta IROC-OROC offset
790 fstringG+="Lx++"; // common slope
791 fstringG+="Lx*isin++"; // delta slope
792 fstringG+="Ly++"; // common slope
793 fstringG+="Ly*isin++"; // delta slope
794 TVectorD vecG[2];
795 TString * strFitG=0;
796 TString * strFitLX=0;
797 //
798 strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideA"+cutAll, chi2,npoints,vecG[0],covar,-1,0, 10000000, kFALSE);
799 chain->SetAlias("tfitGA",strFitG->Data());
800 strFitG->Tokenize("++")->Print();
801 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
802 //
803 strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideC"+cutAll, chi2,npoints,vecG[1],covar,-1,0, 10000000, kFALSE);
804 chain->SetAlias("tfitGC",strFitG->Data());
805 strFitG->Tokenize("++")->Print();
806 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
807 //
808 AliTPCCalPad *padFitG =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[0],vecG[1]);
809 AliTPCCalPad *padFitLX=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[0],vecG[1]);
810 // swap a side and c side
811 AliTPCCalPad *padFitGSwap =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[1],vecG[0]);
812 AliTPCCalPad *padFitLXSwap=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[1],vecG[0]);
813 padFitG->SetName("CEG");
814 padFitLX->SetName("CELX");
815 padFitGSwap->SetName("CEGS");
816 padFitLXSwap->SetName("CELXS");
817 preprocesor->AddComponent(padFitG->Clone());
818 preprocesor->AddComponent(padFitLX->Clone());
819 preprocesor->AddComponent(padFitGSwap->Clone());
820 preprocesor->AddComponent(padFitLXSwap->Clone());
821 preprocesor->DumpToFile("cetmean.root"); // add it to the file
822 //
823 // make local fits
824 //
825 f = new TFile("cetmean.root");
826 chain = (TTree*) f->Get("calPads");
827 AliTPCkalmanAlign::MakeAliasCE(chain);
828 TString fstringL=""; // local fit
829 // // 0. delta common
830 fstringL+="isin++"; // 1. delta IROC-OROC offset
831 fstringL+="Lx++"; // 2. common slope
832 fstringL+="Lx*isin++"; // 3. delta slope
833 fstringL+="Ly++"; // 2. common slope
834 fstringL+="Ly*isin++"; // 3. delta slope
835 TVectorD vecL[36];
836 TVectorD dummy(6);
837 AliTPCCalPad *padFitLCE = new AliTPCCalPad("LocalCE","LocalCE");
838 AliTPCCalPad *padFitTmpCE;
839 for (Int_t isec=0; isec<36; isec++){
840 TCut cutSector=Form("(sector%36)==%d",isec);
841 strFitLX = TStatToolkit::FitPlane(chain,"deltaT-CEG.fElements-CELX.fElements", fstringL.Data(),cutSector+cutAll+"abs(deltaT-CEG.fElements-CELX.fElements)<0.4", chi2,npoints,vecL[isec],covar,-1,0, 10000000, kFALSE);
842 printf("sec=%d\tchi2=%f\n",isec,TMath::Sqrt(chi2/npoints));
843 //
844 TString fitL=Form("((sector%36)==%d)++((sector%36)==%d)*(sector<36)++((sector%36)==%d)*(lx-133)/100.++((sector%36)==%d)*(sector<36)*(lx-133)/100.++((sector%36)==%d)*(ly)/100.++((sector%36)==%d)*(sector<36)*(ly)/100.",isec,isec,isec,isec);
845 if (isec<18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),vecL[isec],dummy);
846 if (isec>=18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),dummy,vecL[isec]);
847 padFitLCE->Add(padFitTmpCE);
848 }
849 //
850 padFitLCE->SetName("CELocal");
851 preprocesor->AddComponent(padFitLCE->Clone());
852 preprocesor->DumpToFile("cetmean.root"); // add it to the file
853 //
854 // write data to array
855 //
856 fFitCEGlobal = new TObjArray(6);
857 fFitCELocal = new TObjArray(6);
858 for (Int_t ipar=0; ipar<6;ipar++){
859 fFitCEGlobal->AddAt(new TVectorD(36),ipar);
860 fFitCELocal->AddAt(new TVectorD(36),ipar);
861 //
862 TVectorD &fvecG = *((TVectorD*)fFitCEGlobal->At(ipar));
863 TVectorD &fvecL = *((TVectorD*)fFitCELocal->At(ipar));
864 //
865 for (Int_t isec=0; isec<36;isec++){
866 fvecL[isec]=vecL[isec][ipar];
867 if (ipar>0){
868 if (isec<18) fvecG[isec]=vecG[0][ipar+2];
869 if (isec>=18) fvecG[isec]=vecG[1][ipar+2];
870 }
871 }
872 }
873 //
874 //
875 //
876}
125d3a38 877
4486a91f 878void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){
879 //
880 // make a aliases of pad variables
881 //
882 chain->SetAlias("side","(-1+(sector%36<18)*2)");
883 chain->SetAlias("sideA","(sector%36<18)");
884 chain->SetAlias("sideC","(sector%36>=18)");
885 chain->SetAlias("isin","(sector<36)");
886 chain->SetAlias("deltaT","CETmean.fElements-PulserTmean.fElements");
887 chain->SetAlias("timeP","PulserTmean.fElements");
888 chain->SetAlias("Gy","(gy.fElements/500.)");
889 chain->SetAlias("Gx","(gx.fElements/500.)");
890 chain->SetAlias("Lx","(lx.fElements-133)/100."); // lx in meters
891 chain->SetAlias("Ly","(ly.fElements)/100.");
892 chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)");
893 chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)");
125d3a38 894}