]>
Commit | Line | Data |
---|---|---|
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 | ||
84 | TObject *palign=0; | |
85 | TTree * treeP=0; | |
86 | TTree * treeM=0; | |
87 | TTree * tree0=0; | |
88 | TTree * treePZ=0; | |
89 | TTree * treeMZ=0; | |
90 | TTree * tree0Z=0; | |
91 | AliTPCROC *proc = AliTPCROC::Instance(); | |
92 | AliTPCExBEffectiveSector * geffCorr=0; | |
93 | ||
94 | Double_t GetCorr(Double_t sector, Double_t localX, Double_t kZ,Int_t type); | |
95 | ||
96 | void MakeFits(); | |
97 | void InitTPCalign(); | |
98 | ||
99 | ||
100 | void 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 | ||
116 | void 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 | ||
134 | void 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 | ||
142 | void 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.="<<¶m<< | |
215 | "\n"; | |
216 | } | |
217 | } | |
218 | } | |
219 | } | |
220 | delete pcstream; | |
221 | } | |
222 | ||
223 | ||
224 | ||
225 | void 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 | ||
263 | Double_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 | ||
288 | void 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 | ||
313 | void 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 | ||
387 | void 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 | ||
414 | void 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 | ||
494 | void 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 |