]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/MakeAlignCalPad.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update variable description file
[u/mrichter/AliRoot.git] / TPC / CalibMacros / MakeAlignCalPad.C
CommitLineData
5b7417d2 1/*
2 marian.ivanov@cern.ch
3 Macro to create alignment/distortion maps
4 As a input output of AliTPCcalibAlign is used.
5
6 Algorithm:
7
8 In the setup without the magnetic field the tracks are fitted using the linear track model.
9 ( in y-rphi coordinate the primary vertex is also used as a constrain)
10 Residuals (deltas0 between the track and clusters in Y and in z direction are filled in the 4 D histogram:
11 Delta: phi(180 bins): localX(53 bins): tan(phi): tan(theta)(10 bins)
12
13 Distortion map are extracted form the residual histograms as a mean value at each bin.
14 Linear fits are then performed for each pad - delta as function of theta
15 Delta Ymeas = offsetY+slopeY*tan(theta)
16 Delta Zmeas = offsetZ+slopeZ*tan(theta)
17
18 Resulting residuals exported into the OCDB are:
19 DeltaY = offsetY
20 DeltaZ = offsetZ
21 DeltaR = slopeZ;
22
23 Example usage:
24
25 make calpad+ make report ps file:
26 aliroot -b -q ~/NimStyle.C ../ConfigCalibTrain.C\(119037\) $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C\(1\)
27 make only report ps file:
28 aliroot -b -q ~/NimStyle.C $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C\(3\)
29 Making fit - iterative procedure - see below:
30
31 gROOT->Macro("~/rootlogon.C");
32 gSystem->Load("libANALYSIS");
33 gSystem->Load("libSTAT");
34 gSystem->Load("libTPCcalib");
35 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros -I$ALICE_ROOT/TPC/TPC -I$ALICE_ROOT/STAT");
36 .L $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C+
37 //load OCDB
38 gROOT->Macro("../ConfigCalibTrain.C(119037)");
39 //2
40 InitTPCalign();
41 MakeFits(); // this is logn proceure 30 minutes
55afc657 42
43 //UpdateOCDB(0,AliCDBRunRange::Infinity());
5b7417d2 44 //
45 gROOT->Macro("~/NimStyle.C")
46
47*/
48#if !defined(__CINT__) || defined(__MAKECINT__)
49#include "TH1D.h"
50#include "TH2F.h"
51#include "THnSparse.h"
52#include "TVectorD.h"
53#include "TTreeStream.h"
54#include "TFile.h"
55#include "AliTPCcalibAlign.h"
56#include "AliTPCROC.h"
57#include "AliTPCCalROC.h"
58#include "AliTPCCalPad.h"
59#include "TF1.h"
60#include "TH2.h"
61#include "TH3.h"
62#include "TROOT.h"
63#include "TProfile.h"
64#include "AliTPCPreprocessorOnline.h"
65#include "AliTPCcalibDB.h"
66#include "AliTPCkalmanAlign.h"
67#include "TPostScript.h"
68#include "TLegend.h"
69#include "TCut.h"
70#include "TCanvas.h"
71#include "TStyle.h"
72#include "AliTPCExBEffectiveSector.h"
73#include "AliTPCComposedCorrection.h"
74//
75#include "AliCDBMetaData.h"
76#include "AliCDBId.h"
77#include "AliCDBManager.h"
78#include "AliCDBStorage.h"
79#include "AliCDBEntry.h"
80#include "TStatToolkit.h"
81#endif
82
83
84TObject *palign=0;
85TTree * treeP=0;
86TTree * treeM=0;
87TTree * tree0=0;
88TTree * treePZ=0;
89TTree * treeMZ=0;
90TTree * tree0Z=0;
91AliTPCROC *proc = AliTPCROC::Instance();
92AliTPCExBEffectiveSector * geffCorr=0;
93
94Double_t GetCorr(Double_t sector, Double_t localX, Double_t kZ,Int_t type);
95
96void MakeFits();
97void InitTPCalign();
98
99
100void MakeAlignCalPad(Int_t mode){
101 //
102 // Make AlignCalpad and make report
103 //
104 gSystem->Load("libANALYSIS");
105 gSystem->Load("libSTAT");
106 gSystem->Load("libTPCcalib");
107 printf("Make report mode\t%d\n", mode);
108 if (mode<1) {
109 InitTPCalign();
110 MakeFits();
111 }
112 printf("Make report\n");
113}
114
115
116void InitTPCalign(){
117 //
118 // read the TPC alignment
119 //
120 TFile fcalib("CalibObjects.root");
121 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
122 if (array){
123 palign = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
124 }else{
125 palign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
126 }
127 if (!palign){
128 TFile fcalib2("TPCAlignObjects.root");
129 palign = ( AliTPCcalibAlign *)fcalib2.Get("alignTPC");
130 }
131}
132
133
134void MakeFits(){
135 AliTPCcalibAlign * align = (AliTPCcalibAlign *)palign;
136 THnSparse * hdY = align->GetClusterDelta(0);
137 THnSparse * hdZ = align->GetClusterDelta(1);
138 AliTPCExBEffectiveSector::MakeResidualMap(hdY,"clusterDY.root");
139 AliTPCExBEffectiveSector::MakeResidualMap(hdZ,"clusterDZ.root");
140}
141
142void FitdY(TTree * tree){
143 //
144 //
145 //
146
147 AliTPCROC * roc = AliTPCROC::Instance();
148 Double_t xquadrant = roc->GetPadRowRadii(36,53);
149 Double_t xIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
150
151 //
152 tree->SetAlias("diffS","(sector-int(sector)-0.5)");
153 tree->SetAlias("diffQ",Form("(localX-%f)",xquadrant));
154 tree->SetAlias("diffIO",Form("(localX-%f)",xIO));
155 tree->SetAlias("iroc",Form("(localX<%f)",xIO)); // bool IROC
156 tree->SetAlias("dqLR0",Form("(localX<%f)*(-1+(diffS<0)*2)",xIO)); //sign LeftRight
157 tree->SetAlias("dqLR2",Form("(localX>%f)*(-1+(diffS<0)*2)",xquadrant)); //sign LeftRight up
158 //
159 tree->SetAlias("diffIFC","abs(localX-80)");
160 tree->SetAlias("diffOFC","abs(250-localX)");
161 tree->SetAlias("drift","(1-abs(kZ*localX)/250.)");
162 //
163 tree->SetAlias("dIFC2",Form("drift/(1+abs(diffIFC^2/(%f^2)))",10.));
164 tree->SetAlias("dIFC",Form("drift/(1+abs(diffIFC/(%f)))",5.));
165 tree->SetAlias("dOFC2",Form("drift/(1+abs(diffOFC^2/(%f^2)))",10.));
166 tree->SetAlias("dOFC",Form("drift/(1+abs(diffOFC/(%f)))",5.));
167
168
169
170 TStatToolkit toolkit;
171 Double_t chi2=0;
172 Int_t npoints=0;
173 TVectorD param;
174 TMatrixD covar;
175 TString fstringG="";
176 fstringG+="(iroc)++"; //1 iroc shift
177 fstringG+="(dqLR0)++";
178 fstringG+="(dqLR2)++";
179 fstringG+="(diffIO*iroc)++"; // iroc rotation
180 fstringG+="(diffIO*dqLR0)++";
181 fstringG+="(diffQ*dqLR2)++";
182 //
183 fstringG+="(dIFC+dIFC2)++"; //9
184 fstringG+="(diffS*(dIFC+dIFC2))++";
185 fstringG+="(diffS^2*(dIFC+dIFC2))++";
186 fstringG+="(dIFC-dIFC2)++"; //9
187 fstringG+="(diffS*(dIFC-dIFC2))++";
188 fstringG+="(diffS^2*(dIFC-dIFC2))++";
189 //
190 //
191 fstringG+="(dOFC+dOFC2)++"; //9
192 fstringG+="(diffS*(dOFC+dOFC2))++";
193 fstringG+="(diffS^2*(dOFC+dOFC2))++";
194 fstringG+="(dOFC-dOFC2)++"; //9
195 fstringG+="(diffS*(dOFC-dOFC2))++";
196 fstringG+="(diffS^2*(dOFC-dOFC2))++";
197 //
198 //
199 //
200 TTreeSRedirector *pcstream=new TTreeSRedirector("dyFit.root");
201 {for (Int_t isec=0; isec<18; isec++){
202 {for (Int_t iside=0; iside<=1; iside++){
203 TCut cutS=Form("abs(sector-%d.5)<0.5&&kZ*%d>0&&entries>2000",isec,(iside>0)?-1:1);
204 TString *strFitG = TStatToolkit::FitPlane(tree,"mean", fstringG.Data(),cutS, chi2,npoints,param,covar,-1,0, 10000000, kTRUE);
205 tree->SetAlias("fitG",strFitG->Data());
206 tree->Draw("mean-fitG",cutS,"");
207 //
208 printf("%d\t%d\t%f\n",iside,isec,TMath::Sqrt(chi2/npoints));
209 (*pcstream)<<"fitDy"<<
210 "iside="<<iside<<
211 "sec="<<isec<<
212 "chi2="<<chi2<<
213 "npoints="<<npoints<<
214 "p.="<<&param<<
215 "\n";
216 }
217 }
218 }
219 }
220 delete pcstream;
221}
222
223
224
225void DumpDerivative(TH3 * his){
226 //
227 //
228 //
229 Int_t nbins0=his->GetXaxis()->GetNbins();
230 Int_t nbins1=his->GetYaxis()->GetNbins();
231 Int_t nbins2=his->GetZaxis()->GetNbins();
232 TTreeSRedirector *pcstream = new TTreeSRedirector("aaa.root");
233 for (Int_t ibin0=1; ibin0<nbins0; ibin0++)
234 for (Int_t ibin1=1; ibin1<nbins1; ibin1++)
235 for (Int_t ibin2=1; ibin2<nbins2; ibin2++){
236 Float_t x= his->GetXaxis()->GetBinCenter(ibin0);
237 Float_t y= his->GetYaxis()->GetBinCenter(ibin1);
238 Float_t z= his->GetZaxis()->GetBinCenter(ibin2);
239 Float_t c= his->GetBinContent(ibin0,ibin1,ibin2);
240 Float_t c0M=his->GetBinContent(ibin0-1,ibin1,ibin2);
241 Float_t c0P=his->GetBinContent(ibin0+1,ibin1,ibin2);
242 Float_t c1M=his->GetBinContent(ibin0,ibin1-1,ibin2);
243 Float_t c1P=his->GetBinContent(ibin0,ibin1+1,ibin2);
244 Float_t c2M=his->GetBinContent(ibin0,ibin1,ibin2-1);
245 Float_t c2P=his->GetBinContent(ibin0,ibin1,ibin2+1);
246 printf("%f\t%f\t%f\t%f\n",x,y,z,c);
247 (*pcstream)<<"aaa"<<
248 "x="<<x<<
249 "y="<<y<<
250 "z="<<z<<
251 "c="<<c<<
252 "c0M="<<c0M<<
253 "c0P="<<c0P<<
254 "c1M="<<c1M<<
255 "c1P="<<c1P<<
256 "c2M="<<c2M<<
257 "c2P="<<c2P<<
258 "\n";
259 }
260 delete pcstream;
261}
262
263Double_t GetCorr(Double_t sector, Double_t localX, Double_t kZ, Int_t type){
264 //
265 // calculate the correction at given position - check the geffCorr
266 //
267 Double_t phi=sector*TMath::Pi()/9.;
268 Double_t gx = localX*TMath::Cos(phi);
269 Double_t gy = localX*TMath::Sin(phi);
270 Double_t gz = localX*kZ;
271 Int_t nsector=(gz>0) ? 0:18;
272 //
273 //
274 //
275 Float_t distPoint[3]={gx,gy,gz};
276 geffCorr->DistortPoint(distPoint, nsector);
277 Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
278 Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
279 Double_t phi0=TMath::ATan2(gy,gx);
280 Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
281 if (type==0) return r1-r0;
282 if (type==1) return (phi1-phi0)*r0;
283 return phi1-phi0;
284}
285
286
287
288void LoadDistortionTrees(){
289 //
290 // Load distortion tree
291 //
292 TFile *fp = new TFile("clusterDYPlus.root");
293 TFile *fm = new TFile("clusterDYMinus.root");
294 TFile *f0 = new TFile("clusterDY0.root");
295 TFile *fpz= new TFile("clusterDZPlus.root");
296 TFile *fmz= new TFile("clusterDZMinus.root");
297 TFile *f0z=new TFile("clusterDZ0.root");
298 treeP=(TTree*)fp->Get("delta");
299 treeM=(TTree*)fm->Get("delta");
300 tree0=(TTree*)f0->Get("delta");
301 treePZ=(TTree*)fpz->Get("delta");
302 treeMZ=(TTree*)fmz->Get("delta");
303 tree0Z=(TTree*)f0z->Get("delta");
304 tree0->AddFriend(treeP,"P");
305 tree0->AddFriend(treeM,"M");
306 tree0->AddFriend(treePZ,"PZ");
307 tree0->AddFriend(treeMZ,"MZ");
308 tree0->AddFriend(tree0Z,"Z");
309 tree0->SetMarkerStyle(25);
310 tree0->SetMarkerSize(0.4);
311}
312
313void UpdateEffSectorOCDB(){
314 //
315 // Incremeantal update ot the correction maps
316 // corrections on top of previous corrections
317 //
318 TFile fp("clusterDYPlus.root");
319 TFile fm("clusterDYMinus.root");
320 TFile f0("clusterDY0.root");
321 TFile f0z("clusterDZ0.root");
322
323 TH3F *his3D0=(TH3F*)f0.Get("his3D");
324 TH3F *his3DP=(TH3F*)fp.Get("his3D");
325 TH3F *his3DM=(TH3F*)fm.Get("his3D");
326 TH3F *his3DZ=(TH3F*)f0z.Get("his3D");
327 TH3F *hisR=0; //half of difference (Plus-Minus) scaled by c1
328 hisR=(TH3F*)his3DP->Clone();
329 hisR->Add(his3DM,-1);
330 hisR->Scale(0.5);
331 hisR->Scale(1/0.3); // c1 factor to be added there //is sign correct?
332 //
333 //
334 AliTPCExBEffectiveSector *effSector = new AliTPCExBEffectiveSector;
335 effSector->SetName("EffSector");
336 effSector->SetTitle("EffSector");
337 effSector->fCorrectionRPhi=(TH3F*)his3D0->Clone();
338 effSector->fCorrectionR=(TH3F*)hisR->Clone();
339 effSector->fCorrectionZ=(TH3F*)his3DZ->Clone();
340 geffCorr=effSector;
55afc657 341 AliTPCCorrection::AddVisualCorrection(geffCorr,0);
5b7417d2 342 //
343 //
344 //
345 AliTPCcalibDB * calib = AliTPCcalibDB::Instance();
346 TObjArray *arr = calib->GetTPCComposedCorrectionArray();
347 {for (Int_t i=0; i<arr->GetEntries(); i++){
348 AliTPCComposedCorrection * corr=(AliTPCComposedCorrection*)arr->At(i);
349 if (!corr) continue;
350 AliTPCExBEffectiveSector * corrSec=( AliTPCExBEffectiveSector *)corr->GetCorrections()->FindObject("EffSector");
351 if (i==0) geffCorr=(AliTPCExBEffectiveSector *)corrSec;
352 if (!corrSec) corr->GetCorrections()->Add(effSector->Clone());
353 if (corrSec){
354 if (corrSec->fCorrectionRPhi) corrSec->fCorrectionRPhi->Add(effSector->fCorrectionRPhi);
355 else corrSec->fCorrectionRPhi=(TH3F*)effSector->fCorrectionRPhi->Clone();
356 if (corrSec->fCorrectionR) corrSec->fCorrectionR->Add(effSector->fCorrectionR);
357 else corrSec->fCorrectionR=(TH3F*)effSector->fCorrectionR->Clone();
358 if (corrSec->fCorrectionZ) corrSec->fCorrectionR->Add(effSector->fCorrectionZ);
359 else corrSec->fCorrectionZ=(TH3F*)effSector->fCorrectionZ->Clone();
360 }
55afc657 361 if (i==0) AliTPCCorrection::AddVisualCorrection(corrSec,1);
5b7417d2 362 }}
363// make OCDB entries
364 TString userName=gSystem->GetFromPipe("echo $USER");
365 TString date=gSystem->GetFromPipe("date");
366 TString ocdbStorage="local:////";
367 ocdbStorage+=gSystem->GetFromPipe("pwd");
368 ocdbStorage+="/OCDB";
369 //
370 AliCDBMetaData *metaData= new AliCDBMetaData();
371 metaData->SetObjectClassName("TObjArray");
372 metaData->SetResponsible("Marian Ivanov");
373 metaData->SetBeamPeriod(1);
374 metaData->SetAliRootVersion("05-27-04"); //root version
375 metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
376 AliCDBId* id1=NULL;
377 id1=new AliCDBId("TPC/Calib/Correction", 0, AliCDBRunRange::Infinity());
378 AliCDBStorage* gStorage = 0;
379
380 gStorage=AliCDBManager::Instance()->GetStorage((ocdbStorage+"Update").Data());
381 gStorage->Put(arr, (*id1), metaData);
382
383
384}
385
386
387void DrawDiff(){
388 TFile f0("../mergeField0/mean.root");
389 TFile fP("../mergePlus/mean.root");
390 TFile fM("../mergeMinus/mean.root");
391 //
392 TTree *itsdy0=(TTree*)f0.Get("ITSdy");
393 TTree *itsdyP=(TTree*)fP.Get("ITSdy");
394 TTree *itsdyM=(TTree*)fM.Get("ITSdy");
395 TTree *vdy0=(TTree*)f0.Get("Vertexdy");
396 TTree *vdyP=(TTree*)fP.Get("Vertexdy");
397 TTree *vdyM=(TTree*)fM.Get("Vertexdy");
398 itsdy0->SetMarkerStyle(25);
399 itsdy0->SetMarkerSize(0.3);
400 itsdy0->AddFriend(itsdyP,"P");
401 itsdy0->AddFriend(itsdyM,"M");
402 itsdy0->AddFriend(vdy0,"V");
403 itsdy0->AddFriend(vdyP,"VP");
404 itsdy0->AddFriend(vdyM,"VM");
405 itsdy0->SetMarkerStyle(25);
406 itsdy0->SetMarkerSize(0.4);
407
408}
409
410
411
412
413
414void MakePlotDeltaZ(){
415 //
416 //
417 //
418 TCut cut="entries>500&&PZ.entries>500&&MZ.entries>500";
419 TCanvas *cA = new TCanvas("deltaZA","deltaZA",900,700);
420 TCanvas *cC = new TCanvas("deltaZC","deltaZC",900,700);
421 TCanvas *cARef = new TCanvas("deltaZARef","deltaZARef",1000,800);
422 TCanvas *cCRef = new TCanvas("deltaZCRef","deltaZCRef",1000,800);
423 //
424 TH1::AddDirectory(0);
425 TH1 * his=0;
426 cA->Divide(4,5);
427 for (Int_t isec=0; isec<18; isec++){
428 cA->cd(isec+1);
429 TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
430 tree0->Draw("Z.mean*10.:localX:abs(kZ)",cutS+"entries>4000&&P.entries>200&&localX<250&&kZ>0&&abs(kZ)<1","colz");
431 his=(TH1*)tree0->GetHistogram()->Clone();
432 his->GetXaxis()->SetTitle("r (cm)");
433 his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
434 his->GetZaxis()->SetTitle("tan($theta)");
435 his->Draw("colz");
436 }
437 cA->cd(19);
438
439 cC->Divide(4,5);
440 for (Int_t isec=0; isec<18; isec++){
441 cC->cd(isec+1);
442 TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
443 tree0->Draw("Z.mean*10.:localX:abs(kZ)",cutS+"entries>4000&&P.entries>200&&localX<250&&kZ<0&&abs(kZ)<1","colz");
444 his=(TH1*)tree0->GetHistogram()->Clone();
445 his->GetXaxis()->SetTitle("r (cm)");
446 his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
447 his->GetZaxis()->SetTitle("tan($theta)");
448 his->Draw("colz");
449
450 }
451
452 cARef->Divide(4,5);
453 for (Int_t isec=0; isec<18; isec++){
454 cARef->cd(isec+1);
455 TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
456 tree0->Draw("(PZ.mean-MZ.mean)*10.:localX:abs(kZ)",cutS+cut+"kZ>0&&abs(kZ)<1","colz");
457 his=(TH1*)tree0->GetHistogram()->Clone();
458 his->GetXaxis()->SetTitle("r (cm)");
459 his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
460 his->GetZaxis()->SetTitle("tan($theta)");
461 his->Draw("colz");
462
463 }
464
465 cCRef->Divide(4,5);
466 for (Int_t isec=0; isec<18; isec++){
467 cCRef->cd(isec+1);
468 TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
55afc657 469 tree0->Draw("(PZ.mean-MZ.mean)*10.:localX:abs(kZ)",cutS+cut+"kZ<0&&abs(kZ)<1","colz");
5b7417d2 470 his=(TH1*)tree0->GetHistogram()->Clone();
471 his->GetXaxis()->SetTitle("r (cm)");
472 his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
473 his->GetZaxis()->SetTitle("tan($theta)");
474 his->Draw("colz");
475 }
476
477 TPostScript *ps = new TPostScript("distortionZ.ps", 112);
478 ps->NewPage();
479 cA->Draw();
480 ps->NewPage();
481 cA->Draw();
482 ps->NewPage();
483 cC->Draw();
484 ps->NewPage();
485 cARef->Draw();
486 ps->NewPage();
487 cCRef->Draw();
488 ps->Close();
489}
490
491
492
493
494void MakeAlign(){
55afc657 495 //
496 // make sector alignment - using Kalman filter method -AliTPCkalmanAlign
497 // Combined information is used, mean residuals are minimized:
498 //
499 // 1. TPC-TPC sector alignment
500 // 2. TPC-ITS alignment
501 // 3. TPC vertex alignment
502 //
5b7417d2 503 TFile fcalib("../mergeField0/TPCAlignObjects.root");
504 AliTPCcalibAlign * align = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
505 TFile f0("../mergeField0/mean.root");
5b7417d2 506
507 //
55afc657 508 TTree *itsdy=(TTree*)f0.Get("ITSdy");
509 TTree *itsdp=(TTree*)f0.Get("ITSdsnp");
510 TTree *itsdz=(TTree*)f0.Get("ITSdz");
511 TTree *itsdt=(TTree*)f0.Get("ITSdtheta");
512 //
513 TTree *vdy=(TTree*)f0.Get("Vertexdy");
514 TTree *vds=(TTree*)f0.Get("Vertexdsnp");
515 TTree *vdz=(TTree*)f0.Get("Vertexdz");
516 TTree *vdt=(TTree*)f0.Get("Vertexdtheta");
5b7417d2 517
55afc657 518 itsdy->AddFriend(itsdp,"Snp");
519 itsdy->AddFriend(itsdz,"Z");
520 itsdy->AddFriend(itsdt,"T");
521 //
522 itsdy->AddFriend(vdy,"V");
523 itsdy->AddFriend(vds,"VSnp");
524 itsdy->AddFriend(vdz,"VZ");
525 itsdy->AddFriend(vdt,"VT");
526 itsdy->SetMarkerStyle(25);
527 itsdy->SetMarkerSize(0.4);
5b7417d2 528
529 TCut cutQ="entries>500&&abs(theta)<0.8&&abs(snp)<0.2";
530 TH1F his1("hdeltaY1","hdeltaY1",100,-0.5,0.5);
531 TMatrixD vecAlign(72,1);
532 TMatrixD covAlign(72,72);
533 AliTPCkalmanAlign::BookAlign1D(vecAlign,covAlign,0,0.0005);
534 TVectorD vecITSY(72);
535 TVectorD vecITSS(72);
536 TVectorD vecVS(72);
55afc657 537 TVectorD vecITSTan(72);
538 TVectorD vecVTan(72);
5b7417d2 539 {for (Int_t isec0=0; isec0<36; isec0++){
540 Double_t phi0=(isec0%18+0.5)*TMath::Pi()/9.;
541 if (phi0>TMath::Pi()) phi0-=TMath::TwoPi();
542 Int_t iside0=(isec0%36<18)? 0:1;
543 TCut cutSector=Form("abs(%f-phi)<0.14",phi0);
544 TCut cutSide = (iside0==0)? "theta>0":"theta<0";
55afc657 545 itsdy->Draw("mean",cutQ+cutSector+cutSide);
546 Double_t meanITSY=itsdy->GetHistogram()->GetMean()/83.6;
5b7417d2 547 vecITSY[isec0]=meanITSY;
548 vecITSY[isec0+36]=meanITSY;
55afc657 549 itsdy->Draw("Snp.mean",cutQ+cutSector+cutSide);
550 Double_t meanITSS=itsdy->GetHistogram()->GetMean();
5b7417d2 551 vecITSS[isec0]=meanITSS;
552 vecITSS[isec0+36]=meanITSS;
55afc657 553 itsdy->Draw("VSnp.mean",cutQ+cutSector+cutSide);
554 Double_t meanVS=itsdy->GetHistogram()->GetMean();
5b7417d2 555 vecVS[isec0]=meanVS;
556 vecVS[isec0+36]=meanVS;
557 }
558 }
559 AliTPCROC * roc = AliTPCROC::Instance();
560 Double_t fXIROC = (roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(0,roc->GetNRows(0)-1))*0.5;
561 Double_t fXOROC = (roc->GetPadRowRadii(36,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
562 Double_t fXmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
563 Double_t fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
564
565 TTreeSRedirector *pcstream=new TTreeSRedirector("combAlign.root");
566
567 {
568 for (Int_t iter=0; iter<2; iter++){
569 for (Int_t isec0=0; isec0<72; isec0++){
570 for (Int_t isec1=0; isec1<72; isec1++){
571 TH1 * his = align->GetHisto(AliTPCcalibAlign::kY,isec0,isec1);
572 TH1 * hisPhi = align->GetHisto(AliTPCcalibAlign::kPhi,isec0,isec1);
573 if (!his) continue;
574 if (his->GetEntries()<100) continue;
575 Double_t xref=fXIO;
576 if (isec0<36&&isec1<36) xref=fXIROC;
577 if (isec0>=36&&isec1>=36) xref=fXOROC;
578 Double_t meanTPC=his->GetMean()/xref;
579 Double_t meanTPCPhi=hisPhi->GetMean();
580 Double_t meanITS0=vecITSY[isec0];
581 Double_t meanITS1=vecITSY[isec1];
582 Double_t meanITSS0=vecITSS[isec0];
583 Double_t meanITSS1=vecITSS[isec1];
584 Double_t meanVS0=vecVS[isec0];
585 Double_t meanVS1=vecVS[isec1];
586 AliTPCkalmanAlign::UpdateAlign1D(meanTPC, 0.001, isec0,isec1, vecAlign,covAlign);
587 AliTPCkalmanAlign::UpdateAlign1D(meanTPCPhi, 0.001, isec0,isec1, vecAlign,covAlign);
588 AliTPCkalmanAlign::UpdateAlign1D(meanITS1-meanITS0, 0.001, isec0,isec1, vecAlign,covAlign);
589 AliTPCkalmanAlign::UpdateAlign1D(meanITSS1-meanITSS0, 0.001, isec0,isec1, vecAlign,covAlign);
590 AliTPCkalmanAlign::UpdateAlign1D(meanVS1-meanVS0, 0.001, isec0,isec1, vecAlign,covAlign);
591 //printf("isec0\t%d\tisec1\t%d\t%f\t%f\t%f\n",isec0,isec1, meanTPC, meanITS0,meanITS1);
592 Double_t kalman0= vecAlign(isec0,0);
593 Double_t kalman1= vecAlign(isec1,0);
594 if (iter>0) (*pcstream)<<"align"<<
595 "iter="<<iter<<
596 "xref="<<xref<<
597 "isec0="<<isec0<<
598 "isec1="<<isec1<<
599 "mTPC="<<meanTPC<<
600 "mTPCPhi="<<meanTPCPhi<<
601 "mITS0="<<meanITS0<<
602 "mITS1="<<meanITS1<<
603 "mITSS0="<<meanITSS0<<
604 "mITSS1="<<meanITSS1<<
605 "mVS0="<<meanVS0<<
606 "mVS1="<<meanVS1<<
607 "k0="<<kalman0<<
608 "k1="<<kalman1<<
609 "\n";
610 }
611 }
612 }
613 }
55afc657 614 pcstream->GetFile()->cd();
615 vecAlign.Write("alignPhiMean");
616 covAlign.Write("alingPhiCovar");
5b7417d2 617 delete pcstream;
618 TFile f("combAlign.root");
619 TTree * treeA = (TTree*)f.Get("align");
620 treeA->SetMarkerStyle(25);
621 treeA->SetMarkerSize(0.5);
622}
623
624
625