]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/FitAlignCombined.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
During simulation: fill STU region w/ non null time sums
[u/mrichter/AliRoot.git] / TPC / CalibMacros / FitAlignCombined.C
CommitLineData
4e093f35 1/*
2 marian.ivanov@cern.ch
3 Macro to create alignment/distortion maps
4 As a input output of AliTPCcalibAlign and AliTPCcalibTime is used.
5 distortion lookup tables are used.
6
7 Input file mean.root with distortion tree expected to be in directory:
8 ../mergeField0/mean.root
9
10
11 The ouput file fitAlignCombined.root contains:
12 1. Resulting (residual) AliTPCCalibMisalignment
13 2. QA fit plots
14
15 Functions documented inside:
16
17 RegisterAlignFunction();
18 MakeAlignFunctionGlobal();
19 MakeAlignFunctionSector();
20
21*/
22
23
24/*
25 Example usage:
26 //
27 .x ~/NimStyle.C
28 gROOT->Macro("~/rootlogon.C");
29 gSystem->Load("libANALYSIS");
30 gSystem->Load("libSTAT");
31 gSystem->Load("libTPCcalib");
32 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros -I$ALICE_ROOT/TPC/TPC -I$ALICE_ROOT/STAT");
33 .L $ALICE_ROOT/TPC/CalibMacros/FitAlignCombined.C+
34 .x ConfigCalibTrain.C(119047)
35 //
36 //
37 FitAlignCombinedCorr();
38
39
40*/
41
42
43
44#if !defined(__CINT__) || defined(__MAKECINT__)
45#include "TH1D.h"
46#include "TH2F.h"
47#include "THnSparse.h"
48#include "TVectorD.h"
49#include "TTreeStream.h"
50#include "TFile.h"
51#include "TChain.h"
52#include "AliTPCcalibAlign.h"
53#include "AliTPCROC.h"
54#include "AliTPCCalROC.h"
55#include "AliTPCCalPad.h"
56#include "TF1.h"
57#include "TH2.h"
58#include "TH3.h"
59#include "TROOT.h"
60#include "TProfile.h"
61#include "AliTPCPreprocessorOnline.h"
62#include "AliTPCcalibDB.h"
63#include "AliTPCkalmanAlign.h"
64#include "TPostScript.h"
65#include "TLegend.h"
66#include "TCut.h"
67#include "TCanvas.h"
68#include "TStyle.h"
69#include "AliLog.h"
70#include "AliTPCExBEffectiveSector.h"
71#include "AliTPCComposedCorrection.h"
72#include "AliTPCCalibGlobalMisalignment.h"
73//
74#include "AliCDBMetaData.h"
75#include "AliCDBId.h"
76#include "AliCDBManager.h"
77#include "AliCDBStorage.h"
78#include "AliCDBEntry.h"
79#include "TStatToolkit.h"
80#include "AliAlignObjParams.h"
81#include "AliTPCParam.h"
82#endif
83
84
85void DrawFitQA();
86
87AliTPCROC *proc = AliTPCROC::Instance();
88AliTPCCalibGlobalMisalignment *combAlignOCDBOld=0; // old misalignment from OCDB - used for reconstruction
89AliTPCCalibGlobalMisalignment *combAlignOCDBNew=0; // new misalignment
90AliTPCCalibGlobalMisalignment *combAlignGlobal=0; // delta misalignment - global part
91AliTPCCalibGlobalMisalignment *combAlignLocal=0; // delta misalignment - sector/local part
92//
93
94TTreeSRedirector *pcstream= 0; // Workspace to store the fit parameters and QA histograms
95//
96TChain * chain=0;
97TChain * chainZ=0; // trees with z measurement
98TTree * tree=0;
99TTree *itsdy=0;
100TTree *itsdyP=0;
101TTree *itsdyM=0;
102TTree *itsdp=0;
103TTree *itsdphiP=0;
104TTree *itsdphiM=0;
105
106TTree *itsdz=0;
107TTree *itsdt=0;
108TTree *tofdy=0;
109TTree *trddy=0;
110//
111TTree *vdy=0;
112TTree *vdyP=0;
113TTree *vdyM=0;
114TTree *vds=0;
115TTree *vdz=0;
116TTree *vdt=0;
117
118TCut cutS="entries>1000&&abs(snp)<0.2&&abs(theta)<1.";
119
120
121void RegisterAlignFunction(){
122 //
123 // Register primitive alignment components.
124 // Linear conbination of primitev forulas used for fit
125 // The nominal delta 1 mm in shift and 1 mrad in rotation
126 // Primitive formulas registeren in AliTPCCoreection::AddvisualCorrection
127 // 0 - deltaX
128 // 1 - deltaY
129 // 2 - deltaZ
130 // 3 - rot0 (phi)
131 // 4 - rot1 (theta)
132 // 5 - rot2
133 TGeoHMatrix matrixX;
134 TGeoHMatrix matrixY;
135 TGeoHMatrix matrixZ;
136 TGeoRotation rot0;
137 TGeoRotation rot1;
138 TGeoRotation rot2; //transformation matrices
139 TGeoRotation rot90; //transformation matrices
140 matrixX.SetDx(0.1); matrixY.SetDy(0.1); matrixZ.SetDz(0.1); //1 mm translation
141 rot0.SetAngles(0.001*TMath::RadToDeg(),0,0);
142 rot1.SetAngles(0,0.001*TMath::RadToDeg(),0);
143 rot2.SetAngles(0,0,0.001*TMath::RadToDeg());
144 //how to get rot02 ?
145 rot90.SetAngles(0,90,0);
146 rot2.MultiplyBy(&rot90,kTRUE);
147 rot90.SetAngles(0,-90,0);
148 rot2.MultiplyBy(&rot90,kFALSE);
149 AliTPCCalibGlobalMisalignment *alignRot0 =new AliTPCCalibGlobalMisalignment;
150 alignRot0->SetAlignGlobal(&rot0);
151 AliTPCCalibGlobalMisalignment *alignRot1=new AliTPCCalibGlobalMisalignment;
152 alignRot1->SetAlignGlobal(&rot1);
153 AliTPCCalibGlobalMisalignment *alignRot2=new AliTPCCalibGlobalMisalignment;
154 alignRot2->SetAlignGlobal(&rot2);
155 //
156 AliTPCCalibGlobalMisalignment *alignTrans0 =new AliTPCCalibGlobalMisalignment;
157 alignTrans0->SetAlignGlobal(&matrixX);
158 AliTPCCalibGlobalMisalignment *alignTrans1=new AliTPCCalibGlobalMisalignment;
159 alignTrans1->SetAlignGlobal(&matrixY);
160 AliTPCCalibGlobalMisalignment *alignTrans2=new AliTPCCalibGlobalMisalignment;
161 alignTrans2->SetAlignGlobal(&matrixZ);
162 AliTPCCorrection::AddVisualCorrection(alignTrans0 ,0);
163 AliTPCCorrection::AddVisualCorrection(alignTrans1 ,1);
164 AliTPCCorrection::AddVisualCorrection(alignTrans2 ,2);
165 AliTPCCorrection::AddVisualCorrection(alignRot0 ,3);
166 AliTPCCorrection::AddVisualCorrection(alignRot1 ,4);
167 AliTPCCorrection::AddVisualCorrection(alignRot2 ,5);
168 TObjArray arrAlign(6);
169 arrAlign.AddAt(alignTrans0->Clone(),0);
170 arrAlign.AddAt(alignTrans1->Clone(),1);
171 arrAlign.AddAt(alignTrans2->Clone(),2);
172 arrAlign.AddAt(alignRot0->Clone(),3);
173 arrAlign.AddAt(alignRot1->Clone(),4);
174 arrAlign.AddAt(alignRot2->Clone(),5);
175 //combAlign.SetCorrections((TObjArray*)arrAlign.Clone());
176}
177
178AliTPCCalibGlobalMisalignment * MakeAlignFunctionGlobal(TVectorD paramYGlobal){
179 //
180 // Take a fit parameters and make a combined correction
181 // 1. Take the common part
182 // 3. Make combined AliTPCCalibGlobalMisalignment - register it
183 // 4. Compare the aliases with fit values - IT is OK
184 //
185 AliTPCCalibGlobalMisalignment *alignGlobal =new AliTPCCalibGlobalMisalignment;
186 TGeoHMatrix matGlobal; // global parameters
187 TGeoHMatrix matDelta; // delta A side - C side
188 //
189 TGeoHMatrix matrixX;
190 TGeoHMatrix matrixY;
191 TGeoHMatrix matrixZ;
192 TGeoRotation rot0;
193 TGeoRotation rot1;
194 TGeoRotation rot2; //transformation matrices
195 TGeoRotation rot90; //transformation matrices
196 //
197 //
198 matrixX.SetDx(0.1*paramYGlobal[1]);
199 matrixY.SetDy(0.1*paramYGlobal[2]);
200 rot0.SetAngles(0.001*TMath::RadToDeg()*paramYGlobal[3],0,0);
201 rot1.SetAngles(0,0.001*TMath::RadToDeg()*paramYGlobal[4],0);
202 rot2.SetAngles(0,0,0.001*TMath::RadToDeg()*paramYGlobal[5]);
203 rot90.SetAngles(0,90,0);
204 rot2.MultiplyBy(&rot90,kTRUE);
205 rot90.SetAngles(0,-90,0);
206 rot2.MultiplyBy(&rot90,kFALSE);
207 matGlobal.Multiply(&matrixX);
208 matGlobal.Multiply(&matrixY);
209 matGlobal.Multiply(&rot0);
210 matGlobal.Multiply(&rot1);
211 matGlobal.Multiply(&rot2);
212 alignGlobal->SetAlignGlobal((TGeoMatrix*)matGlobal.Clone());
213 AliTPCCorrection::AddVisualCorrection(alignGlobal ,100);
214 //
215 /*
216 chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame");
217 chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame");
218
219 chain->Draw("deltaG:fitYGlobal",cutS+"type==0");
220 chain->Draw("deltaG:fitYGlobal",cutS+"type==2");
221
222 */
223
224 return alignGlobal;
225}
226
227
228AliTPCCalibGlobalMisalignment * MakeAlignFunctionSector(TVectorD paramYLocal){
229 //
230 // Take a fit parameters and make a combined correction:
231 // Only delta local Y and delta phi are fitted - not sensityvity for other parameters
232 // Algorithm:
233 // 1. Loop over sectors
234 // 2. Make combined AliTPCCalibGlobalMisalignment - register it
235 // 3. Compare the aliases with fit values - IT is OK
236 //
237 AliTPCCalibGlobalMisalignment *alignLocal =new AliTPCCalibGlobalMisalignment;
238 TGeoHMatrix matrixX;
239 TGeoHMatrix matrixY;
240 TGeoHMatrix matrixZ;
241 TGeoRotation rot0;
242 TGeoRotation rot1;
243 TGeoRotation rot2; //transformation matrices
244 TGeoRotation rot90; //transformation matrices
245 TObjArray * array = new TObjArray(72);
246 TVectorD vecLY(72);
247 TVectorD vecGY(72);
248 TVectorD vecGX(72);
249 TVectorD vecPhi(72);
250 TVectorD vecSec(72);
251 //
252 //
253 {for (Int_t isec=0; isec<72; isec++){
254 Int_t sector=isec%18;
255 Int_t offset=sector*4+1;
256 if ((isec%36)>=18) offset+=2;
257 Double_t phi= (Double_t(sector)+0.5)*TMath::Pi()/9.;
258 //
259 Double_t dly = paramYLocal[offset];
260 Double_t dphi= paramYLocal[offset+1];
261 Double_t dgx = TMath::Cos(phi)*0+TMath::Sin(phi)*dly;
262 Double_t dgy = TMath::Sin(phi)*0-TMath::Cos(phi)*dly;
263 vecSec[isec]=isec;
264 vecLY[isec]=dly;
265 vecGX[isec]=dgx;
266 vecGY[isec]=dgy;
267 vecPhi[isec]=dphi;
268 //
269 matrixX.SetDx(dgx);
270 matrixY.SetDy(dgy);
271 rot0.SetAngles(0.001*TMath::RadToDeg()*dphi,0,0);
272 TGeoHMatrix matrixSector; // global parameters
273 matrixSector.Multiply(&matrixX);
274 matrixSector.Multiply(&matrixY);
275 matrixSector.Multiply(&rot0);
276 array->AddAt(matrixSector.Clone(),isec);
277 }}
278 alignLocal->SetAlignSectors(array);
279 AliTPCCorrection::AddVisualCorrection(alignLocal,101);
280 //
281 //
282 TGraph * graphLY = new TGraph(72,vecSec.GetMatrixArray(), vecLY.GetMatrixArray());
283 TGraph * graphGX = new TGraph(72,vecSec.GetMatrixArray(), vecGX.GetMatrixArray());
284 TGraph * graphGY = new TGraph(72,vecSec.GetMatrixArray(), vecGY.GetMatrixArray());
285 TGraph * graphPhi = new TGraph(72,vecSec.GetMatrixArray(), vecPhi.GetMatrixArray());
286 graphLY->SetMarkerStyle(25); graphGX->SetMarkerStyle(25); graphGY->SetMarkerStyle(25); graphPhi->SetMarkerStyle(25);
287 graphLY->SetName("DeltaLY"); graphLY->SetTitle("#Delta_{ly} (mm)");
288 graphGX->SetName("DeltaGX"); graphGX->SetTitle("#Delta_{gx} (mm)");
289 graphGY->SetName("DeltaGY"); graphGY->SetTitle("#Delta_{gy} (mm)");
290 graphPhi->SetName("DeltaPhi"); graphPhi->SetTitle("#Delta_{phi} (mrad)");
291 //
292 graphLY->Write("grDeltaLY");
293 graphGX->Write("grDeltaGX");
294 graphGY->Write("grDeltaGY");
295 graphPhi->Write("grDeltaPhi");
296 // Check:
297 /*
298 graphLY.Draw("alp")
299 chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame");
300 chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame");
301
302 graphPhi.Draw("alp")
303 chain->Draw("1000*(mean-deltaG):int(sector)",cutS+"type==2&&refX==0&&theta>0","profsame");
304 chain->Draw("1000*(mean-deltaG):int(sector+18)",cutS+"type==2&&refX==0&&theta<0","profsame");
305
306 //
307 chain->Draw("deltaL:fitYLocal",cutS+"type==2&&refX==0",""); // phi shift - OK
308 chain->Draw("deltaL:fitYLocal",cutS+"type==0&&refX==0",""); // OK
309 chain->Draw("deltaL:fitYLocal",cutS+"type==0",""); // OK
310
311 */
312
313 return alignLocal;
314}
315
316
317
318
319
320void LoadTrees(){
321 //
322 // make sector alignment - using Kalman filter method -AliTPCkalmanAlign
323 // Combined information is used, mean residuals are minimized:
324 //
325 // 1. TPC-ITS alignment
326 // 2. TPC vertex alignment
327 //
328 TFile *f0 = new TFile("../mergeField0/mean.root");
329 TFile *fP= new TFile("../mergePlus/mean.root");
330 TFile *fM= new TFile("../mergeMinus/mean.root");
331 //
332 itsdy=(TTree*)f0->Get("ITSdy");
333 itsdyP=(TTree*)fP->Get("ITSdy");
334 itsdyM=(TTree*)fM->Get("ITSdy");
335 itsdp=(TTree*)f0->Get("ITSdsnp");
336 itsdphiP=(TTree*)fP->Get("ITSdsnp");
337 itsdphiM=(TTree*)fM->Get("ITSdsnp");
338
339 itsdz=(TTree*)f0->Get("ITSdz");
340 itsdt=(TTree*)f0->Get("ITSdtheta");
341 tofdy=(TTree*)f0->Get("TOFdy");
342 trddy=(TTree*)f0->Get("TRDdy");
343 //
344 vdy=(TTree*)f0->Get("Vertexdy");
345 vdyP=(TTree*)fP->Get("Vertexdy");
346 vdyM=(TTree*)fM->Get("Vertexdy");
347 vds=(TTree*)f0->Get("Vertexdsnp");
348 vdz=(TTree*)f0->Get("Vertexdz");
349 vdt=(TTree*)f0->Get("Vertexdtheta");
350 chain = new TChain("mean","mean");
351 chainZ = new TChain("mean","mean");
352 chain->AddFile("../mergeField0/mean.root",10000000,"ITSdy");
353 chain->AddFile("../mergeField0/mean.root",10000000,"ITSdsnp");
354 chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdy");
355 chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdsnp");
356 chain->AddFile("../mergeField0/mean.root",10000000,"TOFdy");
357 chain->AddFile("../mergeField0/mean.root",10000000,"TRDdy");
358 //
359 chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdz");
360 chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdtheta");
361 chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdz");
362 chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdtheta");
363
364 //
365 itsdy->AddFriend(itsdp,"Phi");
366 itsdy->AddFriend(itsdphiP,"PhiP");
367 itsdy->AddFriend(itsdphiM,"PhiM");
368 itsdy->AddFriend(itsdz,"Z");
369 itsdy->AddFriend(itsdt,"T");
370 itsdy->AddFriend(itsdyP,"YP");
371 itsdy->AddFriend(itsdyM,"YM");
372 //
373 itsdy->AddFriend(vdy,"V");
374 itsdy->AddFriend(vdyP,"VP");
375 itsdy->AddFriend(vdyP,"VM");
376 itsdy->AddFriend(tofdy,"TOF");
377 itsdy->AddFriend(trddy,"TRD");
378 itsdy->AddFriend(vds,"VPhi");
379 itsdy->AddFriend(vdz,"VZ");
380 itsdy->AddFriend(vdt,"VT");
381 itsdy->AddFriend(tofdy,"TOF.");
382 itsdy->SetMarkerStyle(25);
383 itsdy->SetMarkerSize(0.4);
384 tree=itsdy;
385 tree->SetAlias("side","(-1+2*(theta>0))");
386 chain->SetAlias("side","(-1+2*(theta>0))");
387 chain->SetMarkerStyle(25);
388 chain->SetMarkerSize(0.5);
389}
390
391
392void FitAlignCombinedCorr(){
393 //
394 // Fit Global X and globalY shift at vertex and at ITS
395 //
396 RegisterAlignFunction();
397 LoadTrees();
398 combAlignOCDBOld = AliTPCCalibGlobalMisalignment::CreateOCDBAlign();
399 //
400
401 pcstream= new TTreeSRedirector("fitAlignCombined.root");
402
403 TStatToolkit toolkit;
404 Double_t chi2=0;
405 Int_t npoints=0;
406 // Fit vectors
407 // global fits
408 TVectorD vecSec(37);
409 TVectorD paramYGlobal;
410 TVectorD paramYLocal;
411 TMatrixD covar;
412 // make Aliases
413 {for (Int_t isec=0; isec<18; isec++){ // sectors
414 tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
415 chain->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
416 }}
417 for (Int_t isec=-1;isec<36; isec++){
418 vecSec[isec+1]=isec;
419 }
420 chain->SetAlias("err","((type==0)*0.1+(type==2)*0.001)");
421 // delta phi
422 chain->SetAlias("dpgx","(0+0)");
423 chain->SetAlias("dpgy","(0+0)");
424 chain->SetAlias("dpgr0","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,3)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,3))/160)");
425 chain->SetAlias("dpgr1","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,4)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,4))/160)");
426 chain->SetAlias("dpgr2","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,5)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,5))/160)");
427 chain->SetAlias("deltaPG","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,100)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,100))/160)");
428 chain->SetAlias("deltaPL","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,101)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,101))/160)");
429 // delta y at 85 cm
430 chain->SetAlias("dygxT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,0)+0)");
431 chain->SetAlias("dygyT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,1)+0)");
432 chain->SetAlias("dyr0T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,3)+0)");
433 chain->SetAlias("dyr1T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,4)+0)");
434 chain->SetAlias("dyr2T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,5)+0)");
435 chain->SetAlias("deltaYGT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,100)+0)");
436 chain->SetAlias("deltaYLT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,101)+0)");
437 //
438 // delta y at reference X
439 chain->SetAlias("dygx","(dygxT+0)"); // due global x shift
440 chain->SetAlias("dygy","(dygyT+0)"); // due global y shift
441 chain->SetAlias("dyr0","(dyr0T+dpgr0*(refX-85.))"); // due rotation 0
442 chain->SetAlias("dyr1","(dyr1T+dpgr1*(refX-85.))"); // due rotation 1
443 chain->SetAlias("dyr2","(dyr2T+dpgr2*(refX-85.))"); // due rotation 2
444 chain->SetAlias("deltaYG","(deltaYGT+deltaPG*(refX-85.))"); // due global distortion
445 chain->SetAlias("deltaYL","(deltaYLT+deltaPL*(refX-85.))"); // due local distortion
446
447 chain->SetAlias("deltaG","(type==0)*deltaYG+(type==2)*deltaPG"); //alias to global function
448 chain->SetAlias("deltaL","(type==0)*deltaYL+(type==2)*deltaPL"); //alias to local function
449 //
450 TString fstringGlobal="";
451 // global part
452 fstringGlobal+="((type==0)*dygx+0)++";
453 fstringGlobal+="((type==0)*dygy+0)++";
454 fstringGlobal+="((type==0)*dyr0+(type==2)*dpgr0)++";
455 fstringGlobal+="((type==0)*dyr1+(type==2)*dpgr1)++";
456 fstringGlobal+="((type==0)*dyr2+(type==2)*dpgr2)++";
457 //
458 // Make global fits
459 //
460 TString *strFitYGlobal = TStatToolkit::FitPlane(chain,"mean:err", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYGlobal,covar,-1,0, 10000000, kTRUE);
461 combAlignGlobal =MakeAlignFunctionGlobal(paramYGlobal);
462 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
463 chain->SetAlias("fitYGlobal",strFitYGlobal->Data());
464 paramYGlobal.Print();
465 paramYGlobal.Write("paramYGlobal");
466 //
467 //
468 //
469 //
470 {for (Int_t isec=0; isec<18; isec++){
471 //A side
472 chain->SetAlias(Form("glyASec%d",isec),Form("((type==0)*sec%d*(theta>0))",isec)); // ly shift
473 chain->SetAlias(Form("gxASec%d",isec),Form("((type==0)*dygx*sec%d*(theta>0))",isec)); // gx shift
474 chain->SetAlias(Form("gyASec%d",isec),Form("((type==0)*dygy*sec%d*(theta>0))",isec)); // gy shift
475 chain->SetAlias(Form("gpASec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta>0)",isec)); // phi rotation
476 //C side
477 chain->SetAlias(Form("glyCSec%d",isec),Form("((type==0)*sec%d*(theta<0))",isec)); // ly shift
478 chain->SetAlias(Form("gxCSec%d",isec),Form("((type==0)*dygx*sec%d*(theta<0))",isec));
479 chain->SetAlias(Form("gyCSec%d",isec),Form("((type==0)*dygy*sec%d*(theta<0))",isec));
480 chain->SetAlias(Form("gpCSec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta<0)",isec)); // phi rotation
481 }}
482
483 TString fstringLocal="";
484 {for (Int_t isec=0; isec<18; isec++){
485 //fstringLocal+=Form("((type==0)*dygx*sec%d*(theta>0))++",isec);
486 // fstringLocal+=Form("(gxASec%d)++",isec);
487 // fstringLocal+=Form("(gyASec%d)++",isec);
488 fstringLocal+=Form("(glyASec%d)++",isec);
489 fstringLocal+=Form("(gpASec%d)++",isec);
490 fstringLocal+=Form("(glyCSec%d)++",isec);
491 //fstringLocal+=Form("(gxCSec%d)++",isec);
492 //fstringLocal+=Form("(gyCSec%d)++",isec);
493 fstringLocal+=Form("(gpCSec%d)++",isec);
494 }}
495
496 TString *strFitYLocal = TStatToolkit::FitPlane(chain,"(mean-deltaG):err", fstringLocal.Data(),cutS+"", chi2,npoints,paramYLocal,covar,-1,0, 10000000, kTRUE);
497 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
498 chain->SetAlias("fitYLocal",strFitYLocal->Data());
499 combAlignLocal =MakeAlignFunctionSector(paramYLocal);
500 //paramYLocal.Print();
501 paramYLocal.Write("paramYLocal");
502 //
503 // scaling - why do we need it?
504 TString fstringScale="";
505 fstringScale+="((type==1)*refX+(type==2))*dsec++";
506 fstringScale+="((type==1)*refX+(type==2))*dsec*cos(phi)++";
507 fstringScale+="((type==1)*refX+(type==2))*dsec*sin(phi)++";
508 fstringScale+="((type==1)*refX+(type==2))*dsec*side++";
509 fstringScale+="((type==1)*refX+(type==2))*dsec*side*cos(phi)++";
510 fstringScale+="((type==1)*refX+(type==2))*dsec*side*sin(phi)++";
511 //
512 //
513 //
514 pcstream->GetFile()->cd();
515 DrawFitQA();
516 //
517 AliTPCCalibGlobalMisalignment * combAlignOCDBNew0 =( AliTPCCalibGlobalMisalignment *)combAlignOCDBOld->Clone();
518 combAlignOCDBNew0->AddAlign(*combAlignGlobal);
519 combAlignOCDBNew0->AddAlign(*combAlignLocal);
520 combAlignOCDBNew= AliTPCCalibGlobalMisalignment::CreateMeanAlign(combAlignOCDBNew0);
521 //
522 combAlignGlobal->SetName("fcombAlignGlobal");
523 combAlignLocal->SetName("fcombAlignLocal");
524 combAlignOCDBOld->SetName("fcombAlignOCDBOld");
525 combAlignOCDBNew->SetName("fcombAlignOCDBNew");
526 combAlignOCDBNew0->SetName("fcombAlignOCDBNew0");
527 //
528 combAlignGlobal->Write("fcombAlignGlobal");
529 combAlignLocal->Write("fcombAlignLocal");
530 combAlignOCDBOld->Write("fcombAlignOCDBOld");
531 combAlignOCDBNew->Write("fcombAlignOCDBNew");
532 combAlignOCDBNew0->Write("fcombAlignOCDBNew0");
533 combAlignOCDBNew->Print();
534 delete pcstream;
535}
536
537void DrawFitQA(){
538 //
539 //
540 //
541 //
542 // MakeQA plot 1D
543 //
544 TCanvas c;
545 c.SetLeftMargin(0.15);
546 chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&refX==0","");
547 chain->GetHistogram()->SetName("TPC_VertexDphiGlobal1D");
548 chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)");
549 chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)");
550 chain->GetHistogram()->Fit("gaus","qnr");
551 chain->GetHistogram()->Write();
552 //
553 chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&abs(refX-85)<2","");
554 chain->GetHistogram()->SetName("TPC_ITSDphiGlobal1D");
555 chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)");
556 chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)");
557 chain->GetHistogram()->Fit("gaus","qnr");
558 chain->GetHistogram()->Write();
559 //
560 chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-85)<2","");
561 chain->GetHistogram()->SetName("TPC_ITSDRPhiGlobal1D");
562 chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)");
563 chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)");
564 chain->GetHistogram()->Fit("gaus","qnr");
565 chain->GetHistogram()->Write();
566 //
567 chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-0)<2","");
568 chain->GetHistogram()->SetName("TPC-VertexDRPhiGlobal1D");
569 chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)");
570 chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)");
571 chain->GetHistogram()->Fit("gaus","qnr");
572 chain->GetHistogram()->Write();
573 //
574 // Make QA plot 3D
575 //
576 chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz");
577 chain->GetHistogram()->SetName("TPC_VertexDphi3DGlobal3D");
578 chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)");
579 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
580 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
581 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
582 chain->GetHistogram()->Draw("colz");
583 chain->GetHistogram()->Write();
584 //
585 chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz");
586 chain->GetHistogram()->SetName("TPC_ITSDphi3DGlobal3D");
587 chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)");
588 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
589 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
590 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
591 chain->GetHistogram()->Draw("colz");
592 chain->GetHistogram()->Write();
593 //
594 chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz");
595 chain->GetHistogram()->SetName("TPC_ITSDRphi3DGlobal3D");
596 chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)");
597 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
598 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
599 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
600 chain->GetHistogram()->Draw("colz");
601 chain->GetHistogram()->Write();
602 //
603 chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz");
604 chain->GetHistogram()->SetName("TPC_VertexDRphi3DGlobal3D");
605 chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)");
606 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
607 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
608 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
609 chain->GetHistogram()->Draw("colz");
610 chain->GetHistogram()->Write();
611 //
612 //
613 //
614 chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz");
615 chain->GetHistogram()->SetName("TPC_VertexDphi3DLocal3D");
616 chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Local Fit)");
617 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
618 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
619 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
620 chain->GetHistogram()->Draw("colz");
621 chain->GetHistogram()->Write();
622 //
623 chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz");
624 chain->GetHistogram()->SetName("TPC_ITSDphi3DLocal3D");
625 chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Local Fit)");
626 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
627 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
628 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
629 chain->GetHistogram()->Draw("colz");
630 chain->GetHistogram()->Write();
631 //
632 chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz");
633 chain->GetHistogram()->SetName("TPC_ITSDRphi3DLocal3D");
634 chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Local Fit)");
635 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
636 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
637 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
638 chain->GetHistogram()->Draw("colz");
639 chain->GetHistogram()->Write();
640 //
641 chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz");
642 chain->GetHistogram()->SetName("TPC_VertexDRphiLocal3D");
643 chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Local Fit)");
644 chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
645 chain->GetHistogram()->GetXaxis()->SetTitle("sector");
646 chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
647 chain->GetHistogram()->Draw("colz");
648 chain->GetHistogram()->Write();
649}
650
651
652
653
654void FitAlignCombined0(){
655 //
656 // Fit Global X and globalY shift at vertex and at ITS
657 //
658
659 TTreeSRedirector *pcstream= new TTreeSRedirector("fitAlignCombined.root");
660
661 TStatToolkit toolkit;
662 Double_t chi2=0;
663 Int_t npoints=0;
664 // Fit vectors
665 // global fits
666 TVectorD vecSec(37);
667 TVectorD paramYVGlobal;
668 TVectorD paramYITSGlobal;
669 TVectorD paramPhiVGlobal;
670 TVectorD paramPhiITSGlobal;
671 // local fits
672 TVectorD paramYVLocal;
673 TVectorD paramPhiVLocal;
674 TVectorD paramYITSLocal;
675 TVectorD paramPhiITSLocal;
676 TMatrixD covar;
677 // make Aliases
678 {for (Int_t isec=0; isec<18; isec++){ // sectors
679 tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
680 }}
681 for (Int_t isec=-1;isec<36; isec++){
682 vecSec[isec+1]=isec;
683 }
684 //
685 TString fstringGlobal="";
686 fstringGlobal+="side++";
687 fstringGlobal+="theta++";
688 fstringGlobal+="cos(phi)*(theta>0)++"; // GX - A side
689 fstringGlobal+="sin(phi)*(theta>0)++"; // GY - A side
690 fstringGlobal+="cos(phi)*(theta<0)++"; // GX - C side
691 fstringGlobal+="sin(phi)*(theta<0)++"; // GY - C side
692 fstringGlobal+="cos(phi)*(theta)++"; // theta - get rid of rotation
693 fstringGlobal+="sin(phi)*(theta)++"; // theta - get rid of rotation
694 //
695 // Local Fit
696 //
697 TString fstringLocal="";
698 {for (Int_t sector=0; sector<18; sector++){
699 fstringLocal+=Form("((sec%d)*theta>0)++",sector);
700 fstringLocal+=Form("((sec%d)*theta<0)++",sector);
701 }}
702 //
703 // Make global fits
704 //
705 TString *strFitYVGlobal = TStatToolkit::FitPlane(tree,"V.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYVGlobal,covar,-1,0, 10000000, kFALSE);
706 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
707 tree->SetAlias("fitYVGlobal",strFitYVGlobal->Data());
708 tree->Draw("V.mean:fitYVGlobal",cutS);
709 //
710 TString *strFitPhiVGlobal = TStatToolkit::FitPlane(tree,"VPhi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiVGlobal,covar,-1,0, 10000000, kFALSE);
711 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
712 tree->SetAlias("fitPhiVGlobal",strFitPhiVGlobal->Data());
713 tree->Draw("VPhi.mean:fitPhiVGlobal",cutS);
714 //
715 TString *strFitYITSGlobal = TStatToolkit::FitPlane(tree,"mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYITSGlobal,covar,-1,0, 10000000, kFALSE);
716 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
717 tree->SetAlias("fitYITSGlobal",strFitYITSGlobal->Data());
718 tree->Draw("mean:fitYITSGlobal",cutS);
719 TString *strFitPhiITSGlobal = TStatToolkit::FitPlane(tree,"Phi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiITSGlobal,covar,-1,0, 10000000, kFALSE);
720 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
721 tree->SetAlias("fitPhiITSGlobal",strFitPhiITSGlobal->Data());
722 tree->Draw("Phi.mean:fitPhiITSGlobal",cutS);
723 //
724 // Residuals to the global fit
725 //
726 tree->SetAlias("dyVG", "V.mean-fitYVGlobal");
727 tree->SetAlias("dphiVG","(VPhi.mean-fitPhiVGlobal)");
728 tree->SetAlias("dyITSG","mean-fitYITSGlobal");
729 tree->SetAlias("dphiITSG","(Phi.mean-fitPhiITSGlobal)");
730
731 TString *strFitYVLocal = TStatToolkit::FitPlane(tree,"dyVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYVLocal,covar,-1,0, 10000000, kFALSE);
732 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
733 tree->SetAlias("fitYVLocal",strFitYVLocal->Data());
734 tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz");
735 //
736 TString *strFitPhiVLocal = TStatToolkit::FitPlane(tree,"dphiVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiVLocal,covar,-1,0, 10000000, kFALSE);
737 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
738 tree->SetAlias("fitPhiVLocal",strFitPhiVLocal->Data());
739 tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz");
740 //
741 //
742 //
743 TString *strFitYITSLocal = TStatToolkit::FitPlane(tree,"dyITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYITSLocal,covar,-1,0, 10000000, kFALSE);
744 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
745 tree->SetAlias("fitYITSLocal",strFitYITSLocal->Data());
746 tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
747 //
748 TString *strFitPhiITSLocal = TStatToolkit::FitPlane(tree,"dphiITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiITSLocal,covar,-1,0, 10000000, kFALSE);
749 printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
750 tree->SetAlias("fitPhiITSLocal",strFitPhiITSLocal->Data());
751 tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
752
753 //
754 TH1D *hisA = 0;
755 TH1D *hisC = 0;
756 TVectorD vSec(36);
757 TVectorD vecDPhi(36);
758 TVectorD vecDLY(36);
759 TVectorD vecDGX(36);
760 TVectorD vecDGY(36);
761
762 //
763 tree->SetAlias("phiMean","(fitPhiVLocal+fitPhiITSLocal+fitPhiVGlobal+fitPhiITSGlobal)*0.5");
764 tree->SetAlias("yMeanV","((fitYITSLocal+fitYITSGlobal-85*phiMean)+(fitYVLocal+fitYVGlobal))*0.5");
765 tree->SetAlias("yMeanITS","((fitYITSLocal+fitYITSGlobal)+(fitYVLocal+fitYVGlobal+85*phiMean))*0.5");
766
767
768 tree->Draw("phiMean:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof");
769 hisA = (TH1D*) tree->GetHistogram()->Clone();
770 tree->Draw("phiMean:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof");
771 hisC = (TH1D*) tree->GetHistogram()->Clone();
772
773 for (Int_t isec=0; isec<36; isec++){
774 vSec[isec]=isec;
775 if (isec<18) vecDPhi[isec]=hisA->GetBinContent(isec%18+1);
776 if (isec>=18) vecDPhi[isec]=hisC->GetBinContent(isec%18+1);
777 }
778 tree->Draw("yMeanV:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof");
779 hisA = (TH1D*) tree->GetHistogram()->Clone();
780 tree->Draw("yMeanV:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof");
781 hisC = (TH1D*) tree->GetHistogram()->Clone();
782 //
783 for (Int_t isec=0; isec<36; isec++){
784 if (isec<18) vecDLY[isec]=hisA->GetBinContent(isec%18+1);
785 if (isec>=18) vecDLY[isec]=hisC->GetBinContent(isec%18+1);
786 vecDGX[isec]=-vecDLY[isec]*TMath::Sin(TMath::Pi()*(isec+0.5)/9.);
787 vecDGY[isec]= vecDLY[isec]*TMath::Cos(TMath::Pi()*(isec+0.5)/9.);
788 }
789
790
791
792
793 //
794 // Store results
795 //
796 {(*pcstream)<<"fitLocal"<<
797 //global fits
798 "sec.="<<&vSec<<
799 "pYVGlobal.="<<&paramYVGlobal<<
800 "pYITSGlobal.="<<&paramYITSGlobal<<
801 "pPhiVGlobal.="<<&paramPhiVGlobal<<
802 "pPhiITSGlobal.="<<&paramPhiITSGlobal<<
803 // local fits
804 "pYVLocal.="<<&paramYVLocal<<
805 "pPhiVLocal.="<<&paramPhiVLocal<<
806 "pYITSLocal.="<<&paramYITSLocal<<
807 "pPhiITSLocal.="<<&paramPhiITSLocal<<
808 //
809 // combined
810 "vDPhi.="<<&vecDPhi<<
811 "vDLY.="<<&vecDLY<<
812 "vDGX.="<<&vecDGX<<
813 "vDGY.="<<&vecDGY<<
814 "\n";
815 }
816 paramYVGlobal.Write("paramYVGlobal");
817 paramYITSGlobal.Write("paramYITSGlobal");
818 paramPhiVGlobal.Write("paramPhiVGlobal");
819 paramPhiITSGlobal.Write("paramPhiITSGlobal");
820 // local fits
821 paramYVLocal.Write("paramYVLocal");
822 paramPhiVLocal.Write("paramPhiVLocal");
823 paramYITSLocal.Write("paramYITSLocal");
824 paramPhiITSLocal.Write("paramPhiITSLocal");
825 vecDPhi.Write("vecDPhi");
826 vecDLY.Write("vecDLY");
827 vecDGX.Write("vecDGX");
828 vecDGY.Write("vecDGY");
829
830
831
832 tree->Draw("dyVG:sector:abs(theta)",cutS+"theta>0","colz");
833 tree->GetHistogram()->Write("fitYVGlobal");
834 tree->Draw("dphiVG:sector:abs(theta)",cutS+"theta>0","colz");
835 tree->GetHistogram()->Write("fitPhiVGlobal");
836 tree->Draw("dyITSG:sector:abs(theta)",cutS+"theta>0","colz");
837 tree->GetHistogram()->Write("fitYITSGlobal");
838 tree->Draw("dphiITSG:sector:abs(theta)",cutS+"theta>0","colz");
839 tree->GetHistogram()->Write("fitPhiITSGlobal");
840 //
841 tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz");
842 tree->GetHistogram()->Write("fitYVLocal");
843 tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz");
844 tree->GetHistogram()->Write("fitPhiVLocal");
845 tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
846 tree->GetHistogram()->Write("fitYITSLocal");
847 tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
848 tree->GetHistogram()->Write("fitPhiITSLocal");
849 delete pcstream;
850}
851
852void FitAlignCombined(){
853 //
854 //
855 // make sector alignment - using Kalman filter method -AliTPCkalmanAlign
856 // Combined information is used, mean residuals are minimized:
857 //
858 // 1. TPC-TPC sector alignment
859 // 2. TPC-ITS alignment
860 // 3. TPC vertex alignment
861 TFile fcalib("../mergeField0/TPCAlignObjects.root");
862 AliTPCcalibAlign * align = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
863
864 TCut cutQ="entries>1000&&abs(theta)<0.5&&abs(snp)<0.2&&YP.entries>50&&YM.entries>50";
865 TH1F his1("hdeltaY1","hdeltaY1",100,-0.5,0.5);
866 TMatrixD vecAlign(72,1);
867 TMatrixD covAlign(72,72);
868 TMatrixD vecAlignY(72,1);
869 TMatrixD covAlignY(72,72);
870 TMatrixD vecAlignTheta(72,1);
871 TMatrixD covAlignTheta(72,72);
872 TMatrixD vecAlignZ(72,1);
873 TMatrixD covAlignZ(72,72);
874 AliTPCkalmanAlign::BookAlign1D(vecAlign,covAlign,0,0.002);
875 AliTPCkalmanAlign::BookAlign1D(vecAlignY,covAlignY,0,0.1);
876 AliTPCkalmanAlign::BookAlign1D(vecAlignTheta,covAlignTheta,0,0.002);
877 AliTPCkalmanAlign::BookAlign1D(vecAlignZ,covAlignZ,0,0.1);
878 TVectorD vecITSY(72);
879 TVectorD vecITSYPM(72);
880 TVectorD vecITSPhi(72);
881 TVectorD vecITSPhiPM(72);
882 TVectorD vecVY(72);
883 TVectorD vecVS(72);
884 TVectorD vecITSTheta(72);
885 TVectorD vecVTheta(72);
886 {for (Int_t isec0=0; isec0<36; isec0++){
887 Double_t phi0=(isec0%18+0.5)*TMath::Pi()/9.;
888 if (phi0>TMath::Pi()) phi0-=TMath::TwoPi();
889 Int_t iside0=(isec0%36<18)? 0:1;
890 TCut cutSector=Form("abs(%f-phi)<0.14",phi0);
891 TCut cutSide = (iside0==0)? "theta>0":"theta<0";
892 //
893 itsdy->Draw("mean",cutQ+cutSector+cutSide);
894 Double_t meanITSY=itsdy->GetHistogram()->GetMean();
895 vecITSY[isec0]=meanITSY;
896 vecITSY[isec0+36]=meanITSY;
897 //
898 itsdy->Draw("(YP.mean+YM.mean)*0.5",cutQ+cutSector+cutSide);
899 Double_t meanITSYPM=itsdy->GetHistogram()->GetMean();
900 vecITSYPM[isec0]=meanITSYPM;
901 vecITSYPM[isec0+36]=meanITSYPM;
902 //
903 itsdy->Draw("Phi.mean",cutQ+cutSector+cutSide);
904 Double_t meanITSPhi=itsdy->GetHistogram()->GetMean();
905 vecITSPhi[isec0]=meanITSPhi;
906 vecITSPhi[isec0+36]=meanITSPhi;
907 //
908 itsdy->Draw("(PhiP.mean+PhiM.mean)*0.5",cutQ+cutSector+cutSide);
909 Double_t meanITSPhiPM=itsdy->GetHistogram()->GetMean();
910 vecITSPhiPM[isec0]=meanITSPhiPM;
911 vecITSPhiPM[isec0+36]=meanITSPhiPM;
912 //
913 //
914 itsdy->Draw("VPhi.mean",cutQ+cutSector+cutSide);
915 Double_t meanVS=itsdy->GetHistogram()->GetMean();
916 vecVS[isec0]=meanVS;
917 vecVS[isec0+36]=meanVS;
918 //
919 itsdy->Draw("V.mean",cutQ+cutSector+cutSide);
920 Double_t meanVY=itsdy->GetHistogram()->GetMean();
921 vecVY[isec0]=meanVY;
922 vecVY[isec0+36]=meanVY;
923 //
924 itsdy->Draw("T.mean",cutQ+cutSector+cutSide);
925 Double_t meanITST=itsdy->GetHistogram()->GetMean();
926 vecITSTheta[isec0]=meanITST;
927 vecITSTheta[isec0+36]=meanITST;
928 //
929 itsdy->Draw("VT.mean",cutQ+cutSector+cutSide);
930 Double_t meanVT=itsdy->GetHistogram()->GetMean();
931 vecVTheta[isec0]=meanVT;
932 vecVTheta[isec0+36]=meanVT;
933 }
934 }
935 AliTPCROC * roc = AliTPCROC::Instance();
936 Double_t fXIROC = (roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(0,roc->GetNRows(0)-1))*0.5;
937 Double_t fXOROC = (roc->GetPadRowRadii(36,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
938 Double_t fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
939
940 TTreeSRedirector *pcstream=new TTreeSRedirector("combAlign.root");
941
942 {
943 for (Int_t iter=0; iter<5; iter++){
944 for (Int_t isec0=0; isec0<72; isec0++){
945 for (Int_t isec1=0; isec1<72; isec1++){
946 //if (isec0%36!=isec0%36) continue;
947 if (iter==0 && isec0%36!=isec1%36) continue;
948 if (iter==1 && isec0%18!=isec1%18) continue;
949 TH1 * his = align->GetHisto(AliTPCcalibAlign::kY,isec0,isec1);
950 TH1 * hisPhi = align->GetHisto(AliTPCcalibAlign::kPhi,isec0,isec1);
951 TH1 * hisTheta = align->GetHisto(AliTPCcalibAlign::kTheta,isec0,isec1);
952 TH1 * hisZ = align->GetHisto(AliTPCcalibAlign::kZ,isec0,isec1);
953 Int_t side0=(isec0%36)<18 ? 1:-1;
954 Int_t side1=(isec1%36)<18 ? 1:-1;
955 Double_t weightTPC= 0.002;
956 if (isec0%18==isec1%18) weightTPC=0.0005;
957 if (isec0%36==isec1%36) weightTPC=0.00005;
958 if (side0!=side1) weightTPC*=2.;
959 Double_t weightTPCT= 0.001;
960 if (isec0%18==isec1%18) weightTPC=0.0005;
961 if (isec0%36==isec1%36) weightTPC=0.00005;
962 if (TMath::Abs(isec0%18-isec1%18)==1) weightTPC = 0.0005;
963 if (TMath::Abs(isec0%36-isec1%36)==1) weightTPC = 0.0001;
964
965 //
966 Double_t meanITS0=vecITSY[isec0];
967 Double_t meanITSPM0=vecITSYPM[isec0];
968 Double_t meanITSPhi0=vecITSPhi[isec0];
969 Double_t meanITSPhiPM0=vecITSPhiPM[isec0];
970 Double_t meanVY0=vecVY[isec0];
971 Double_t meanVPhi0=vecVS[isec0];
972 Double_t meanITSTheta0=vecITSTheta[isec0];
973 Double_t meanVTheta0=vecVTheta[isec0];
974 //
975 Double_t meanITS1=vecITSY[isec1];
976 Double_t meanITSPM1=vecITSYPM[isec1];
977 Double_t meanITSPhi1=vecITSPhi[isec1];
978 Double_t meanITSPhiPM1=vecITSPhiPM[isec1];
979 Double_t meanVY1=vecVY[isec1];
980 Double_t meanVPhi1=vecVS[isec1];
981 Double_t meanITSTheta1=vecITSTheta[isec1];
982 Double_t meanVTheta1=vecVTheta[isec1];
983 //
984 Double_t meanPhi0 = (meanITSPhi0+meanVPhi0+2*meanITS0/83.)/4.;
985 Double_t meanPhi1 = (meanITSPhi1+meanVPhi1+2*meanITS1/83.)/4.;
986 //
987 //
988 if (iter>2 &&isec0==isec1){
989 AliTPCkalmanAlign::UpdateAlign1D(-meanITS0/83, 0.001, isec0, vecAlign,covAlign);
990 AliTPCkalmanAlign::UpdateAlign1D(-meanITSPhi0, 0.001, isec0, vecAlign,covAlign);
991 AliTPCkalmanAlign::UpdateAlign1D(-meanVPhi0, 0.001, isec0, vecAlign,covAlign);
992 //
993 AliTPCkalmanAlign::UpdateAlign1D(-meanPhi0, 0.001, isec0, vecAlign,covAlign);
994 AliTPCkalmanAlign::UpdateAlign1D(-meanVTheta0, 0.001, isec0, vecAlignTheta,covAlignTheta);
995 AliTPCkalmanAlign::UpdateAlign1D(-meanITSTheta0, 0.001, isec0, vecAlignTheta,covAlignTheta);
996 }
997 if (iter>2){
998 AliTPCkalmanAlign::UpdateAlign1D(meanPhi1-meanPhi0, 0.001, isec0,isec1, vecAlign,covAlign);
999 AliTPCkalmanAlign::UpdateAlign1D(meanITSPhiPM1-meanITSPhiPM0, 0.001, isec0,isec1, vecAlign,covAlign);
1000 AliTPCkalmanAlign::UpdateAlign1D((meanITSPM1-meanITSPM0)/83., 0.001, isec0,isec1, vecAlign,covAlign);
1001 AliTPCkalmanAlign::UpdateAlign1D(meanVTheta1-meanVTheta0, 0.001, isec0,isec1, vecAlignTheta,covAlignTheta);
1002 AliTPCkalmanAlign::UpdateAlign1D(meanITSTheta1-meanITSTheta0, 0.001, isec0,isec1, vecAlignTheta,covAlignTheta);
1003 }
1004 //
1005 if (!his) continue;
1006 if (!hisPhi) continue;
1007 if (!hisTheta) continue;
1008 if (his->GetEntries()<100) continue;
1009 Double_t xref=fXIO;
1010 if (isec0<36&&isec1<36) xref=fXIROC;
1011 if (isec0>=36&&isec1>=36) xref=fXOROC;
1012 Double_t meanTPC=his->GetMean();
1013 Double_t meanTPCPhi=hisPhi->GetMean();
1014 Double_t meanTPCTheta=hisTheta->GetMean();
1015 Double_t meanTPCZ=hisZ->GetMean();
1016 //
1017 //
1018 Double_t kalman0= vecAlign(isec0,0);
1019 Double_t kalman1= vecAlign(isec1,0);
1020 Double_t kalmanY0= vecAlignY(isec0,0);
1021 Double_t kalmanY1= vecAlignY(isec1,0);
1022 Double_t kalmanTheta0= vecAlignTheta(isec0,0);
1023 Double_t kalmanTheta1= vecAlignTheta(isec1,0);
1024 Double_t kalmanZ0= vecAlignZ(isec0,0);
1025 Double_t kalmanZ1= vecAlignZ(isec1,0);
1026 //
1027 //
1028 AliTPCkalmanAlign::UpdateAlign1D((meanTPC)/xref, weightTPC, isec0,isec1, vecAlign,covAlign);
1029 AliTPCkalmanAlign::UpdateAlign1D(meanTPCPhi, weightTPC*10,isec0,isec1, vecAlign,covAlign);
1030 //
1031 AliTPCkalmanAlign::UpdateAlign1D(meanTPCTheta, weightTPCT, isec0,isec1, vecAlignTheta,covAlignTheta);
1032 if (side0==side1) AliTPCkalmanAlign::UpdateAlign1D(meanTPCZ, weightTPCT*100., isec0,isec1, vecAlignZ,covAlignZ);
1033 //printf("isec0\t%d\tisec1\t%d\t%f\t%f\t%f\n",isec0,isec1, meanTPC, meanITS0,meanITS1);
1034
1035 if (iter>=0) (*pcstream)<<"align"<<
1036 "iter="<<iter<<
1037 "xref="<<xref<< // reference X
1038 "isec0="<<isec0<< // sector number
1039 "isec1="<<isec1<<
1040 "side0="<<side0<<
1041 "side1="<<side1<<
1042 //TPC
1043 "mTPC="<<meanTPC<< // delta Y / xref
1044 "mTPCPhi="<<meanTPCPhi<< // delta Phi
1045 "mTPCZ="<<meanTPCZ<< // delta Z
1046 "mTPCTheta="<<meanTPCTheta<< // delta Theta
1047 //ITS
1048 "mITS0="<<meanITS0<<
1049 "mITS1="<<meanITS1<<
1050 "mITSPhi0="<<meanITSPhi0<<
1051 "mITSPhi1="<<meanITSPhi1<<
1052 //
1053 "mITSPM0="<<meanITSPM0<<
1054 "mITSPM1="<<meanITSPM1<<
1055 "mITSPhiPM0="<<meanITSPhiPM0<<
1056 "mITSPhiPM1="<<meanITSPhiPM1<<
1057 //
1058 "mITSTheta0="<<meanITSTheta0<<
1059 "mITSTheta1="<<meanITSTheta1<<
1060 //Vertex
1061 "mVY0="<<meanVY0<<
1062 "mVY1="<<meanVY1<<
1063 "mVPhi0="<<meanVPhi0<<
1064 "mVPhi1="<<meanVPhi1<<
1065 "mVTheta0="<<meanVTheta0<<
1066 "mVTheta1="<<meanVTheta1<<
1067 // Vertex+ITS mean
1068 "mPhi0="<<meanPhi0<<
1069 "mPhi1="<<meanPhi1<<
1070 //Kalman
1071 "kY0="<<kalmanY0<<
1072 "kY1="<<kalmanY1<<
1073 "kPhi0="<<kalman0<<
1074 "kPhi1="<<kalman1<<
1075 "kTheta0="<<kalmanTheta0<<
1076 "kTheta1="<<kalmanTheta1<<
1077 "kZ0="<<kalmanZ0<<
1078 "kZ1="<<kalmanZ1<<
1079 "\n";
1080 }
1081 }
1082 }
1083 }
1084 pcstream->GetFile()->cd();
1085 vecAlign.Write("alignPhiMean");
1086 covAlign.Write("alingPhiCovar");
1087 vecAlignTheta.Write("alignThetaMean");
1088 covAlignTheta.Write("alingThetaCovar");
1089 vecAlignZ.Write("alignZMean");
1090 covAlignZ.Write("alingZCovar");
1091 delete pcstream;
1092 TFile f("combAlign.root");
1093 TTree * treeA = (TTree*)f.Get("align");
1094 treeA->SetMarkerStyle(25);
1095 treeA->SetMarkerSize(0.5);
1096}
1097
1098
1099
1100
1101
1102void UpdateOCDBAlign(){
1103 //
1104 // Store resulting OCDB entry
1105 // 0. Setup OCDB to get necccessary old entries - not done here
1106 // 1. Get old OCDB entry
1107 // 2. Get delta alignment
1108 // 3. Add delta alignment
1109 // 4. Store new alignment in
1110
1111 AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
1112 TClonesArray * array = (TClonesArray*)entry->GetObject();
1113 Int_t entries = array->GetEntries();
1114 TClonesArray * arrayNew=(TClonesArray*)array->Clone();
1115 TFile f("combAlign.root");
1116 TMatrixD *matPhi=(TMatrixD*)f.Get("alignPhiMean");
1117 //
1118 //
1119 { for (Int_t i=0;i<entries; i++){
1120 //
1121 //
1122 AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
1123 AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i);
1124 Int_t imod;
1125 AliGeomManager::ELayerID ilayer;
1126 alignP->GetVolUID(ilayer, imod);
1127 if (ilayer==AliGeomManager::kInvalidLayer) continue;
1128 Int_t sector=imod;
1129 if (ilayer==AliGeomManager::kTPC2) sector+=36;
1130 Double_t transOld[3], rotOld[3];
1131 alignP->GetTranslation(transOld); // in cm
1132 alignP->GetAngles(rotOld); // in degrees
1133 printf("%d\t%d\t%d\t",ilayer, imod,sector);
1134 printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*matPhi)(sector,0)));
1135 alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]+(*matPhi)(sector,0)*TMath::RadToDeg());
1136 alignPNew->Print("kokot");
1137 alignP->Print("kokot");
1138 }
1139 }
1140
1141
1142
1143 TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate";
1144 TString userName=gSystem->GetFromPipe("echo $USER");
1145 TString date=gSystem->GetFromPipe("date");
1146
1147 AliCDBMetaData *metaData= new AliCDBMetaData();
1148 metaData->SetObjectClassName("TObjArray");
1149 metaData->SetResponsible("Marian Ivanov");
1150 metaData->SetBeamPeriod(1);
1151 metaData->SetAliRootVersion("05-26-02"); //root version
1152 metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
1153
1154 AliCDBId* id1=NULL;
1155 id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity());
1156 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1157 gStorage->Put(arrayNew, (*id1), metaData);
1158
1159 TFile falign("falign.root","recreate");
1160 arrayNew->Write("new");
1161 array->Write("old");
1162 falign.Close();
1163}
1164
1165
1166
1167void UpdateOCDBAlign0(){
1168 //
1169 // Store resulting OCDB entry
1170 // 0. Setup OCDB to get necccessary old entries - not done here
1171 // 1. Get old OCDB entry
1172 // 2. Get delta alignment
1173 // 3. Add delta alignment
1174 // 4. Store new alignment in
1175
1176 AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
1177 TClonesArray * array = (TClonesArray*)entry->GetObject();
1178 Int_t entries = array->GetEntries();
1179 TClonesArray * arrayNew=(TClonesArray*)array->Clone();
1180 TFile f("fitAlignCombined.root");
1181 TVectorD *vecDPhi = (TVectorD*) f.Get("vecDPhi");
1182 TVectorD *vecDGX = (TVectorD*) f.Get("vecDGX");
1183 TVectorD *vecDGY = (TVectorD*) f.Get("vecDGY");
1184
1185 //
1186 //
1187 { for (Int_t i=0;i<entries; i++){
1188 //
1189 //
1190 AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
1191 AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i);
1192 Int_t imod;
1193 AliGeomManager::ELayerID ilayer;
1194 alignP->GetVolUID(ilayer, imod);
1195 if (ilayer==AliGeomManager::kInvalidLayer) continue;
1196 Int_t sector=imod;
1197 if (ilayer==AliGeomManager::kTPC2) sector+=36;
1198 Double_t transOld[3], rotOld[3];
1199 alignP->GetTranslation(transOld); // in cm
1200 alignP->GetAngles(rotOld); // in degrees
1201 printf("%d\t%d\t%d\t",ilayer, imod,sector);
1202 printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*vecDPhi)[sector%36]));
1203 alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]-(*vecDPhi)[sector%36]*TMath::RadToDeg());
1204 alignPNew->SetTranslation(transOld[0]-(*vecDGX)[sector%36],transOld[1]-(*vecDGY)[sector%36], transOld[2]);
1205 alignPNew->Print("kokot");
1206 alignP->Print("kokot");
1207 }
1208 }
1209
1210
1211
1212 TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate";
1213 TString userName=gSystem->GetFromPipe("echo $USER");
1214 TString date=gSystem->GetFromPipe("date");
1215
1216 AliCDBMetaData *metaData= new AliCDBMetaData();
1217 metaData->SetObjectClassName("TObjArray");
1218 metaData->SetResponsible("Marian Ivanov");
1219 metaData->SetBeamPeriod(1);
1220 metaData->SetAliRootVersion("05-26-02"); //root version
1221 metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
1222
1223 AliCDBId* id1=NULL;
1224 id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity());
1225 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1226 gStorage->Put(arrayNew, (*id1), metaData);
1227
1228 TFile falign("falign.root","recreate");
1229 arrayNew->Write("new");
1230 array->Write("old");
1231 falign.Close();
1232}