]>
Commit | Line | Data |
---|---|---|
59fbd3a7 | 1 | #if !defined(__CINT__) || defined(__MAKECINT__) |
2 | #include <iostream> | |
3 | #include "TROOT.h" | |
4 | #include "TSelector.h" | |
5 | #include "TGeoManager.h" | |
6 | #include "TArrayI.h" | |
7 | #include "TText.h" | |
8 | #include "TCanvas.h" | |
9 | #include "TStyle.h" | |
10 | #include "TFile.h" | |
11 | #include "TTree.h" | |
12 | #include "TNtuple.h" | |
13 | #include "TH1.h" | |
14 | #include "TH2.h" | |
15 | #include "AliTRDgeometry.h" | |
16 | #include "AliLog.h" | |
17 | #include "AliAlignObj.h" | |
18 | #include "AliAlignObjParams.h" | |
19 | #include "AliAlignmentTracks.h" | |
20 | #include "AliTrackFitterStraight.h" | |
21 | #include "AliTrackFitterRieman.h" | |
22 | #include "AliTrackFitterKalman.h" | |
23 | #include "AliTrackResidualsLinear.h" | |
24 | #include "AliTrackResidualsFast.h" | |
25 | #include "AliTrackResidualsChi2.h" | |
26 | #include "AliTRDalignment.h" | |
27 | #endif | |
28 | ||
29 | //#include "/misc/misko/root/peakfit.C" | |
30 | ||
31 | AliAlignmentTracks *alt; | |
32 | AliTrackFitter *trf; | |
33 | AliTrackResiduals *res; | |
34 | TNtuple *ntresi; // residuals after minimization | |
35 | TNtuple *ntalpa; // alignment params | |
36 | TNtuple *nteva; // results of minimization | |
37 | double aml[100][40][6]; // results from align_module_layers [iter][layer][phi,z,r,rphi,rz,rr] | |
38 | enum {kStraight, kRieman, kKalman}; // track fitters | |
39 | enum {kLinear, kFast, kChi2}; // minimizer | |
40 | ||
41 | //============================================================================= | |
42 | TNtuple *makeNtuple(AliTrackResiduals *res, char *name, int flag) { | |
43 | ||
44 | // make and fill ntuple containing residuals so they can be visualized | |
45 | // flag=1 means including transformation by AliTrackResiduals::fAlignObj | |
46 | // use flag 0 to see the status before minimization | |
47 | // use flag 1 to see the status after minimization | |
48 | ||
49 | TNtuple *nt = new TNtuple(name,name,"x0:y0:z0:x1:y1:z1:v:dy",320000); | |
50 | AliTrackPoint p0,p1; // 0 - reference; 1 - actual point | |
51 | ||
52 | for (Int_t itrack = 0; itrack < res->GetNFilledTracks(); itrack++) { | |
53 | AliTrackPointArray *volarray; | |
54 | AliTrackPointArray *traarray; | |
55 | if (!res->GetTrackPointArrays(itrack, volarray, traarray)) continue; | |
56 | for (Int_t ipoint = 0; ipoint < volarray->GetNPoints(); ipoint++) { | |
57 | traarray->GetPoint(p0,ipoint); | |
58 | volarray->GetPoint(p1,ipoint); | |
59 | if (flag) res->GetAlignObj()->Transform(p1); | |
60 | AliTrackPoint q0,q1; // in tracking system | |
61 | Double_t alpha = p1.GetAngle(); | |
62 | q0=p0.Rotate(alpha); | |
63 | q1=p1.Rotate(alpha); | |
64 | UShort_t volid = p1.GetVolumeID(); | |
65 | double x0 = q0.GetX(); | |
66 | double y0 = q0.GetY(); | |
67 | double z0 = q0.GetZ(); | |
68 | double x1 = q1.GetX(); | |
69 | double y1 = q1.GetY(); | |
70 | double z1 = q1.GetZ(); | |
71 | double dy = y1-y0; | |
72 | AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid); | |
73 | dy += 0.035*(z1-z0)*pow(-1.0,layer); // pad tilt | |
74 | nt->Fill(x0,y0,z0,x1,y1,z1,volid,dy); | |
75 | } | |
76 | } | |
77 | return nt; | |
78 | } | |
79 | //============================================================================= | |
80 | void filterTracks(AliTrackResiduals *res, double level) { | |
81 | ||
82 | // find the peak in y-residuals; remove outlier tracks | |
83 | ||
84 | // first, find the peak and prepare the track quality cut | |
85 | ||
86 | char *hinam = "filter_tracks_tmp"; | |
87 | TH1D *hist = new TH1D(hinam,hinam,250,-5,5); | |
88 | TNtuple *nt = makeNtuple(res,"kuku",0); | |
89 | nt->Project(hinam,"dy"); | |
90 | delete nt; | |
91 | ||
92 | if (hist->GetEntries()<100) return; | |
93 | ||
94 | // find peak in the histogram | |
95 | Int_t bpeak=(Int_t) hist->GetMaximumBin(); | |
96 | Int_t nbins=(Int_t) hist->GetXaxis()->GetNbins(); | |
97 | Int_t ypeak=(Int_t) hist->GetBinContent(bpeak); | |
98 | if (bpeak==0 || bpeak==nbins) return; | |
99 | ||
100 | // find the edges of the distribution | |
101 | Int_t bleft,brigt; | |
102 | {for (bleft=bpeak; bleft>0; bleft--) { | |
103 | if (hist->GetBinContent(bleft)<level*ypeak) break; | |
104 | }} | |
105 | {for (brigt=bpeak; brigt<nbins; brigt++) { | |
106 | if (hist->GetBinContent(brigt)<level*ypeak) break; | |
107 | }} | |
108 | bleft=bpeak+(bleft-bpeak); | |
109 | brigt=bpeak+(brigt-bpeak); | |
110 | if (bleft<0) bleft=0; | |
111 | if (brigt>nbins) brigt=nbins; | |
112 | Axis_t xleft=hist->GetBinCenter(bleft); | |
113 | Axis_t xrigt=hist->GetBinCenter(brigt); | |
114 | delete hist; | |
115 | ||
116 | // second, make copy of original tracks | |
117 | ||
118 | Int_t ntracks = res->GetNFilledTracks(); | |
119 | if (ntracks<1) return; | |
120 | AliTrackPointArray **tmpvol = new AliTrackPointArray*[ntracks]; | |
121 | AliTrackPointArray **tmptra = new AliTrackPointArray*[ntracks]; | |
122 | for (Int_t itrack = 0; itrack < ntracks; itrack++){ | |
123 | AliTrackPointArray *volarray; | |
124 | AliTrackPointArray *traarray; | |
125 | if (!res->GetTrackPointArrays(itrack, volarray, traarray)) {printf("cannot access track point array"); return;} | |
126 | tmpvol[itrack] = new AliTrackPointArray(*volarray); | |
127 | tmptra[itrack] = new AliTrackPointArray(*traarray); | |
128 | } | |
129 | ||
130 | // third, fire all employees and hire only the good ones | |
131 | ||
132 | res->SetNTracks(ntracks); | |
133 | {for (Int_t itrack = 0; itrack < ntracks; itrack++) { | |
134 | AliTrackPointArray *volarray = tmpvol[itrack]; | |
135 | AliTrackPointArray *traarray = tmptra[itrack]; | |
136 | int good_track=1; | |
137 | {for (Int_t ipoint = 0; ipoint < volarray->GetNPoints(); ipoint++) { | |
138 | AliTrackPoint p0,p1; | |
139 | traarray->GetPoint(p0,ipoint); | |
140 | volarray->GetPoint(p1,ipoint); | |
141 | // res->GetAlignObj()->Transform(p1); | |
142 | AliTrackPoint q0,q1; // in tracking system | |
143 | Double_t alpha = p1.GetAngle(); | |
144 | q0=p0.Rotate(alpha); | |
145 | q1=p1.Rotate(alpha); | |
146 | double y0 = q0.GetY(); | |
147 | double z0 = q0.GetZ(); | |
148 | double y1 = q1.GetY(); | |
149 | double z1 = q1.GetZ(); | |
150 | UShort_t volid = p1.GetVolumeID(); | |
151 | AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid); | |
152 | double dy = y1-y0-0.035*(z1-z0)*pow(-1.0,layer); | |
153 | if (dy<xleft || dy>xrigt) good_track=0; | |
154 | }} | |
155 | if (good_track) res->AddTrackPointArrays(volarray, traarray); | |
156 | }} | |
157 | } | |
158 | //============================================================================= | |
159 | void drawWithTitles(TTree *tr, char *what, char *xtitle, char *ytitle) { | |
160 | tr->Draw(what); | |
161 | tr->GetHistogram()->GetXaxis()->SetTitle(xtitle); | |
162 | tr->GetHistogram()->GetYaxis()->SetTitle(ytitle); | |
163 | tr->GetHistogram()->DrawCopy(); | |
164 | } | |
165 | //============================================================================= | |
166 | TCanvas *show(TNtuple *nt) { | |
167 | printf("%d track points\n",(int) nt->GetEntries()); | |
168 | gStyle->SetOptFit(1); | |
169 | gStyle->SetOptTitle(0); | |
170 | gStyle->SetTitleOffset(3.5,"y"); | |
171 | gStyle->SetTitleOffset(2.5,"x"); | |
172 | gStyle->SetStatX(0.99); | |
173 | gStyle->SetStatY(0.99); | |
174 | gStyle->SetStatW(0.3); | |
175 | gStyle->SetPadLeftMargin(0.15); | |
176 | gStyle->SetPadRightMargin(0.10); | |
177 | TCanvas *c = new TCanvas(); | |
178 | c->Divide(3,2); | |
179 | c->cd(1); drawWithTitles(nt, "y1-y0:y0>>hi(60,-60,60,40,-2,2)", "y_{track} (cm)", "y-y_{track} (cm)"); | |
180 | c->cd(2); drawWithTitles(nt, "y1-y0:x0>>hi(60,280,400,40,-2,2)", "x_{track} (cm)", "y-y_{track} (cm)"); | |
181 | c->cd(3); drawWithTitles(nt, "y1-y0:z0>>hi(70,-350,350,40,-2,2)", "z_{track} (cm)", "y-y_{track} (cm)"); | |
182 | c->cd(4); drawWithTitles(nt, "y1-y0:z1-z0>>hi(50,-10,10,40,-2,2)", "z-z_{track} (cm)", "y-y_{track} (cm)"); | |
183 | c->cd(5); drawWithTitles(nt, "dy>>dyhist(150,-2,4)", "dy +- 0.035*dz", ""); // same binning as in FilterTracks | |
184 | // double par[10]; | |
185 | // peakfit(nt->GetHistogram(),0.1,0.1,par); | |
186 | nt->GetHistogram()->Fit("gaus"); | |
187 | ((TH1F*) gDirectory->Get("dyhist"))->DrawCopy(); | |
188 | return c; | |
189 | } | |
190 | //============================================================================= | |
191 | void init_alt(const char *trapo, int fitter_flag, int res_flag) { | |
192 | ||
193 | // prepare the AliAlignmentTracks object | |
194 | // read space points from file trapo | |
195 | ||
196 | alt = new AliAlignmentTracks(); | |
197 | alt->SetPointsFilename(trapo); | |
198 | //alt->Misalign("kuku.root","TRDAlignObjs"); | |
199 | ||
200 | if (fitter_flag == kStraight) { | |
201 | trf = new AliTrackFitterStraight(); | |
202 | } else if (fitter_flag == kRieman) { | |
203 | trf = new AliTrackFitterRieman(); | |
204 | ((AliTrackFitterRieman *) trf)->SetMaxDelta(6); | |
205 | } else if (fitter_flag == kKalman) { | |
206 | trf = new AliTrackFitterKalman(); | |
207 | ((AliTrackFitterKalman *) trf)->SetMaxChi2(1000); | |
208 | } | |
209 | trf->SetMinNPoints(30); // if working with clusters | |
210 | trf->SetMinNPoints(4); // if working with tracklets | |
211 | alt->SetTrackFitter(trf); | |
212 | ||
213 | if (res_flag == kLinear) { | |
214 | res = new AliTrackResidualsLinear(); | |
215 | ((AliTrackResidualsLinear *)res)->SetRobust(0.7); | |
216 | } else if (res_flag == kFast) { | |
217 | res = new AliTrackResidualsFast(); | |
218 | } else if (res_flag == kChi2) { | |
219 | res = new AliTrackResidualsChi2(); | |
220 | res->FixParameter(0,0); | |
221 | res->FixParameter(1,0); | |
222 | res->FixParameter(2,0); | |
223 | res->FixParameter(3,0); | |
224 | res->FixParameter(4,0); | |
225 | // only rotation around z allowed | |
226 | } | |
227 | res->SetMinNPoints(10); // if working with clusters | |
228 | res->SetMinNPoints(1); // if working with tracklets | |
229 | alt->SetMinimizer(res); | |
230 | ||
231 | if (nteva) delete nteva; | |
232 | nteva = new TNtuple("eva","eva","chamber:module:layer:iteration:p0:p1:p2:p3:p4:p5:fp0:fp1:fp2"); | |
233 | } | |
234 | //============================================================================= | |
235 | void misa_alt(UShort_t volid, Double_t *mis) { | |
236 | ||
237 | // Misalign volid in alt by (local) pa[6]. | |
238 | // Actually, both fAlignObjs and fMisalignObjs are applied in alt. | |
239 | // The proper way to introduce misalignment would be via fMisalignObjs. | |
240 | // Because of lazyness I am using fAlignObj (they exist already). | |
241 | // As a consequence, one should not expect the resulting transformations | |
242 | // to be the inverse of the introduced misalignment, but to be zero. | |
243 | ||
244 | AliAlignObjParams *al = (AliAlignObjParams *) alt->GetAlignObj(volid); | |
245 | al->SetLocalPars(mis[0],mis[1],mis[2],mis[3],mis[4],mis[5]); | |
246 | } | |
247 | //============================================================================= | |
248 | int align_volume(TArrayI volIds, TArrayI volIdRefs, int iterations, | |
249 | double level, double *pa, int showflag) { | |
250 | ||
251 | // fit tracks in volumes volIdRefs | |
252 | // determine residuals in volIds and save on AliTrackResiduals res | |
253 | // find alignment parameters such as to minimize the residuals | |
254 | // repeat the procedure iteration times | |
255 | // return resulting params in pa | |
256 | // if level>0 then suppress tracks from the tails of the distribution | |
257 | // tails being defined as the places where the distribution falls below | |
258 | // level (relative to the peak height) | |
259 | // showflag = 1 show the residuals before minimization | |
260 | // showflag = 2 show the residuals after minimization | |
261 | ||
262 | // align | |
263 | ||
264 | printf("Aligning volume"); | |
265 | if (volIds.GetSize()>1) printf("s"); | |
266 | for (int i=0; i<volIds.GetSize(); i++) printf(" %d (%s)",volIds.At(i),AliGeomManager::SymName(volIds.At(i))); | |
267 | printf(" to reference volume"); | |
268 | if (volIdRefs.GetSize()>1) printf("s"); | |
269 | for (int i=0; i<volIdRefs.GetSize(); i++) printf(" %d (%s)",volIdRefs.At(i),AliGeomManager::SymName(volIdRefs.At(i))); | |
270 | printf("\n"); | |
271 | ||
272 | int resu = alt->AlignVolumes(&volIds, &volIdRefs, AliGeomManager::kTPC1, AliGeomManager::kTOF,iterations); | |
273 | if (!resu) return 0; | |
274 | ||
275 | // repeat taking only peak | |
276 | ||
277 | if (level) { | |
278 | filterTracks(res,level); | |
279 | // alt->AlignVolumes(&volIds, &volIdRefs, AliGeomManager::kTPC1, AliGeomManager::kTOF,iterations); | |
280 | // this would reload res... just do | |
281 | AliAlignObjParams oldal(*res->GetAlignObj()); // old alignment object | |
282 | res->Minimize(); | |
283 | AliAlignObjParams newal(*res->GetAlignObj()); // new alignment object | |
284 | AliAlignObjParams delta(oldal.Inverse()); | |
285 | delta *= newal; // new/old | |
286 | oldal.GetPars(pa,pa+3); printf("oldal %10.3f\n",pa[1]); | |
287 | newal.GetPars(pa,pa+3); printf("newal %10.3f\n",pa[1]); | |
288 | delta.GetPars(pa,pa+3); printf("delta %10.3f\n",pa[1]); | |
289 | if (alt->GetUpdate()) for (Int_t i = 0; i < volIds.GetSize(); i++) { | |
290 | AliAlignObj *alignObj = alt->GetAlignObj(volIds.At(i)); | |
291 | *alignObj *= delta; | |
292 | } | |
293 | } | |
294 | ||
295 | printf("minimizer: %12d tracks\n",res->GetNFilledTracks()); | |
296 | printf("minimizer: %12d degrees of freedom\n",res->GetNdf()); | |
297 | printf("minimizer: %12.2f chi2\n",res->GetChi2()); | |
298 | printf("minimizer: %12.2f normalized chi2\n",res->GetChi2()/res->GetNdf()); | |
299 | ||
300 | {for (int i=0; i<volIds.GetSize(); i++) { | |
301 | // TGeoHMatrix m; | |
302 | // res->GetAlignObj()->GetMatrix(m); | |
303 | // alt->GetAlignObj(volIds.At(i))->SetMatrix(m); // delta matrix rather than final matrix, for debugging | |
304 | alt->GetAlignObj(volIds.At(i))->GetLocalPars(pa,pa+3); | |
305 | printf("updated alignment object for %d %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n", | |
306 | volIds.At(i),pa[0],pa[1],pa[2],pa[3],pa[4],pa[5]); | |
307 | }} | |
308 | ||
309 | // show result | |
310 | ||
311 | if (showflag) { | |
312 | if (ntresi) delete ntresi; | |
313 | ntresi = makeNtuple(res,"kuku",showflag-1); | |
314 | TCanvas *c = show(ntresi); | |
315 | ||
316 | c->cd(6); | |
317 | TText te; | |
318 | te.SetTextFont(43); | |
319 | te.SetTextColor(1); | |
320 | te.SetTextSizePixels(11); | |
321 | double xpos = 0.1; | |
322 | double ypos = 1.0; | |
323 | te.SetTextAlign(13); | |
324 | char buf[2000]; | |
325 | ||
326 | te.DrawTextNDC(xpos,ypos-=0.04,res->GetName()); | |
327 | te.DrawTextNDC(xpos,ypos-=0.07,"Aligning volumes"); | |
328 | for (int i=0; i<volIds.GetSize(); i++) { | |
329 | sprintf(buf,"%d (%s) ",volIds.At(i),AliGeomManager::SymName(volIds.At(i))); | |
330 | te.DrawTextNDC(xpos+0.05,ypos-=0.04,buf); | |
331 | if (volIds.GetSize()>4 && i==1) { // cut it short | |
332 | te.DrawTextNDC(xpos+0.05,ypos-=0.04,"...etc..."); | |
333 | i=volIds.GetSize()-2; | |
334 | } | |
335 | } | |
336 | te.DrawTextNDC(xpos,ypos-=0.04,"to reference volumes"); | |
337 | for (int i=0; i<volIdRefs.GetSize(); i++) { | |
338 | sprintf(buf,"%d (%s) ",volIdRefs.At(i),AliGeomManager::SymName(volIdRefs.At(i))); | |
339 | te.DrawTextNDC(xpos+0.05,ypos-=0.04,buf); | |
340 | } | |
341 | ||
342 | te.DrawTextNDC(xpos,ypos-=0.07,"Result"); | |
343 | double ypos1 = ypos; | |
344 | te.DrawTextNDC(xpos+0.05,ypos1-=0.04,"shift in phi"); | |
345 | te.DrawTextNDC(xpos+0.05,ypos1-=0.04,"shift in z"); | |
346 | te.DrawTextNDC(xpos+0.05,ypos1-=0.04,"shift in r"); | |
347 | te.DrawTextNDC(xpos+0.05,ypos1-=0.04,"tilt in phi"); | |
348 | te.DrawTextNDC(xpos+0.05,ypos1-=0.04,"tilt in z"); | |
349 | te.DrawTextNDC(xpos+0.05,ypos1-=0.04,"tilt in r"); | |
350 | te.SetTextAlign(33); | |
351 | sprintf(buf,"%.3f",pa[0]); te.DrawTextNDC(xpos+0.45,ypos-=0.04,buf); | |
352 | sprintf(buf,"%.3f",pa[1]); te.DrawTextNDC(xpos+0.45,ypos-=0.04,buf); | |
353 | sprintf(buf,"%.3f",pa[2]); te.DrawTextNDC(xpos+0.45,ypos-=0.04,buf); | |
354 | sprintf(buf,"%.4f",pa[3]); te.DrawTextNDC(xpos+0.45,ypos-=0.04,buf); | |
355 | sprintf(buf,"%.4f",pa[4]); te.DrawTextNDC(xpos+0.45,ypos-=0.04,buf); | |
356 | sprintf(buf,"%.4f",pa[5]); te.DrawTextNDC(xpos+0.45,ypos-=0.04,buf); | |
357 | te.SetTextAlign(13); | |
358 | if (showflag==1) te.DrawTextNDC(xpos,ypos-=0.07,"Residuals shown before alignment"); | |
359 | if (showflag==2) te.DrawTextNDC(xpos,ypos-=0.07,"Residuals shown after alignment"); | |
360 | c->Print("c1.ps"); | |
361 | } | |
362 | ||
363 | return res->GetNFilledTracks(); | |
364 | } | |
365 | //============================================================================= | |
366 | int align_volume(int module, int layer0, int layer1, int reflayer0, int reflayer1, | |
367 | int iterations, double level, double *pa, int showflag) { | |
368 | ||
369 | // align the chambers (module,layer0-layer1) to the layers reflayer0-reflayer1 | |
370 | // (but excluding the layers to be aligned) | |
371 | ||
372 | // prepare volume arrays | |
373 | int n = 0; | |
374 | int temp[1000]={}; | |
375 | int sec = module/5; | |
376 | for (int i=reflayer0; i<=reflayer1; i++) { | |
377 | if (i>=layer0 && i<=layer1) continue; | |
378 | //if (i!=reflayer0 && i!=reflayer1) continue; // take only first and last ref layer | |
379 | if (i>=AliGeomManager::kTPC1 && i<=AliGeomManager::kTPC2) { | |
380 | // both z-halves of TPC | |
381 | temp[n++] = AliGeomManager::LayerToVolUID(i,sec); | |
382 | temp[n++] = AliGeomManager::LayerToVolUID(i,sec+18); | |
383 | } | |
384 | if (i>=AliGeomManager::kTRD1 && i<=AliGeomManager::kTRD6) { | |
385 | temp[n++] = AliGeomManager::LayerToVolUID(i,module); | |
386 | } | |
387 | } | |
388 | TArrayI volIdRefs(n); | |
389 | for (int i=0; i<n; i++) volIdRefs.AddAt(temp[i],i); | |
390 | ||
391 | TArrayI volIds(layer1-layer0+1); | |
392 | ||
393 | for (int i=layer0; i<=layer1; i++) volIds.AddAt(AliGeomManager::LayerToVolUID(i,module),i-layer0); | |
394 | ||
395 | return align_volume(volIds,volIdRefs,iterations,level,pa,showflag); | |
396 | } | |
397 | //============================================================================= | |
398 | int align_sm(int sec, int iterations, double level, double *pa, int showflag) { | |
399 | ||
400 | // align the whole supermodule to the TPC | |
401 | // could be unified with the align_volume above | |
402 | ||
403 | // prepare volume arrays | |
404 | TArrayI volIdRefs(4); | |
405 | volIdRefs.AddAt(AliGeomManager::LayerToVolUID(7,sec),0); | |
406 | volIdRefs.AddAt(AliGeomManager::LayerToVolUID(8,sec),1); | |
407 | volIdRefs.AddAt(AliGeomManager::LayerToVolUID(7,sec+18),2); | |
408 | volIdRefs.AddAt(AliGeomManager::LayerToVolUID(8,sec+18),3); | |
409 | ||
410 | TArrayI volIds(30); | |
411 | int n = 0; | |
412 | for (int layer=9; layer<=14; layer++) for (int module=5*sec; module<5*sec+5; module++) | |
413 | volIds.AddAt(AliGeomManager::LayerToVolUID(layer,module),n++); | |
414 | ||
415 | return align_volume(volIds,volIdRefs,iterations,level,pa,showflag); | |
416 | } | |
417 | //============================================================================= | |
418 | void align_layer_modules(int layer, int reflayer0, int reflayer1, | |
419 | int iterations, double level) { | |
420 | ||
421 | // Loop through modules (TRD chambers) within given layer. | |
422 | // Align each module using other layers. | |
423 | // Store result on an ntuple. | |
424 | ||
425 | double a[6]={}; | |
426 | double b[6]={}; | |
427 | ntalpa = new TNtuple("ntalpa","ntalpa","a0:a1:a2:a3:a4:a5:b0:b1:b2:b3:b4:b5"); | |
428 | Int_t maxmodule = AliGeomManager::LayerSize(layer); | |
429 | {for (int i=0; i<maxmodule; i++) { | |
430 | printf("************* module %d of layer %d *************\n",i,layer); | |
431 | int ntra = align_volume(i,layer,layer,reflayer0,reflayer1,iterations,level,b,0); | |
432 | if (ntra) ntalpa->Fill(a[0],a[1],a[2],a[3],a[4],a[5],b[0],b[1],b[2],b[3],b[4],b[5]); | |
433 | }} | |
434 | } | |
435 | //============================================================================= | |
436 | void align_module_layers(int module, int iterations, double level) { | |
437 | ||
438 | // Loop through the 4 inner layers (TRD chambers) within given stack. | |
439 | // Find the alignment needed to match each chamber to the other 5. | |
440 | // Apply all 4 alignments. | |
441 | // Repeat this procedure iteration times. | |
442 | ||
443 | int flayer = AliGeomManager::kTRD1; | |
444 | int llayer = AliGeomManager::kTRD6; | |
445 | char buf[1000]; | |
446 | double par[10]; | |
447 | TH1F *dy=0; | |
448 | for (int i=0; i<iterations; i++) { | |
449 | // go through the layers once | |
450 | {for (int layer=flayer+1; layer<llayer; layer++) { | |
451 | //if (layer==11) continue; | |
452 | align_volume(module, layer, layer, flayer, llayer, 1, level, aml[i][layer], 0); | |
453 | int sector = module/5; | |
454 | int stack = module%5; | |
455 | int det = AliTRDgeometry::GetDetector(layer-AliGeomManager::kTRD1,stack,sector); | |
456 | nteva->Fill(det,module,layer,i, | |
457 | aml[i][layer][0],aml[i][layer][1],aml[i][layer][2], | |
458 | aml[i][layer][3],aml[i][layer][4],aml[i][layer][5], | |
459 | 0,0,0); | |
460 | UShort_t voluid = AliGeomManager::LayerToVolUID(layer,module); | |
461 | alt->GetAlignObj(voluid)->SetLocalPars(0,0,0,0,0,0); | |
462 | if (ntresi) delete ntresi; | |
463 | ntresi = makeNtuple(res,"kuku",0); | |
464 | sprintf(buf,"dy>>hidy%02d_%d_befo(100,-2,2)",i,layer); | |
465 | ntresi->Draw(buf); | |
466 | sprintf(buf,"hidy%02d_%d_befo",i,layer); | |
467 | dy=(TH1F*) gDirectory->Get(buf); | |
468 | // peakfit(dy,0.1,0.1,par); | |
469 | nteva->Fill(det,module,layer,i-0.2,0,0,0,0,0,0, | |
470 | par[0],par[1],par[2]); | |
471 | ||
472 | if (ntresi) delete ntresi; | |
473 | ntresi = makeNtuple(res,"kuku",1); | |
474 | sprintf(buf,"dy>>hidy%02d_%d_afte(100,-2,2)",i,layer); | |
475 | ntresi->Draw(buf); | |
476 | sprintf(buf,"hidy%02d_%d_afte",i,layer); | |
477 | dy=(TH1F*) gDirectory->Get(buf); | |
478 | // peakfit(dy,0.1,0.1,par); | |
479 | dy->DrawCopy(); | |
480 | nteva->Fill(det,module,layer,i+0.2,0,0,0,0,0,0, | |
481 | par[0],par[1],par[2]); | |
482 | }} | |
483 | // update all 6 alignment objects | |
484 | {for (int layer=flayer+1; layer<llayer; layer++) { | |
485 | UShort_t voluid = AliGeomManager::LayerToVolUID(layer,module); | |
486 | alt->GetAlignObj(voluid)->SetLocalPars( | |
487 | aml[i][layer][0], | |
488 | aml[i][layer][1], | |
489 | aml[i][layer][2], | |
490 | aml[i][layer][3], | |
491 | aml[i][layer][4], | |
492 | aml[i][layer][5]); | |
493 | }} | |
494 | } | |
495 | ||
496 | } | |
497 | //============================================================================= | |
498 | void align_module_layers_plot(int iterations, double ymax) { | |
499 | // show histograms produced by align_module_layers | |
500 | ||
501 | int flayer = AliGeomManager::kTRD1; | |
502 | int llayer = AliGeomManager::kTRD6; | |
503 | int ncol = 8; | |
504 | if (2*iterations>ncol) ncol=2*iterations; | |
505 | TCanvas *c = new TCanvas(); | |
506 | /* | |
507 | char *xdes[]={"befo","afte","befo","afte","befo","afte","befo","afte","befo","afte", | |
508 | "befo","afte","befo","afte","befo","afte","befo","afte","befo","afte"}; | |
509 | char *ydes[]={"plane 5","plane 4","plane 3","plane 2","plane 1","plane 0"}; | |
510 | multi(ncol,6, -2,2,0,ymax, 0,0, 4,5, "dy (cm)", "", xdes, ydes, ""); | |
511 | */ | |
512 | c->Divide(ncol,6); | |
513 | char buf[1000]; | |
514 | for (int layer=flayer; layer<=llayer; layer++) { | |
515 | int col=0; | |
516 | for (int i=0; i<iterations; i++) { | |
517 | int row = AliGeomManager::kTRD6-layer; | |
518 | sprintf(buf,"hidy%02d_%d_befo",i,layer); | |
519 | TH1F *hidy = (TH1F*) gDirectory->Get(buf); | |
520 | if (!hidy) break; | |
521 | hidy->SetMaximum(ymax); | |
522 | //printf("drawing %s %ld in pad %d\n",buf,hidy,row*ncol+i+1); | |
523 | c->cd(row*ncol+col+1); hidy->Draw("same"); col++; | |
524 | sprintf(buf,"hidy%02d_%d_afte",i,layer); | |
525 | hidy = (TH1F*) gDirectory->Get(buf); | |
526 | hidy->SetMaximum(ymax); | |
527 | //printf("drawing %s %ld in pad %d\n",buf,hidy,row*ncol+i+1); | |
528 | c->cd(row*ncol+col+1); hidy->Draw("same"); col++; | |
529 | c->Modified(); | |
530 | c->Update(); | |
531 | } | |
532 | } | |
533 | } | |
534 | //============================================================================= | |
535 | void col(TNtuple *nt, int i) { | |
536 | nt->SetLineColor(i); | |
537 | nt->SetMarkerColor(i); | |
538 | } | |
539 | //============================================================================= | |
540 | void align_module_layers_summary_plot() { | |
541 | ||
542 | // didi(1); | |
543 | // fosi(16); | |
544 | TCanvas *c = new TCanvas("c","c"); | |
545 | c->Divide(2,2); | |
546 | TH2D *du0 = new TH2D("du0","du0",20,-0.5,9.5,100,-0.2,0.2); | |
547 | TH2D *du1 = new TH2D("du1","du1",20,-0.5,9.5,100,-0.2,0.2); | |
548 | TH2D *du2 = new TH2D("du2","du2",20,-0.5,9.5,100,0,0.5); | |
549 | du0->SetXTitle("iteration"); | |
550 | du1->SetXTitle("iteration"); | |
551 | du2->SetXTitle("iteration"); | |
552 | du0->SetYTitle("suggested shift in phi (cm)"); | |
553 | du1->SetYTitle("residuals peak position (cm)"); | |
554 | du2->SetYTitle("residuals peak width (cm)"); | |
555 | ||
556 | c->cd(1); du0->Draw(); | |
557 | col(nteva,1); nteva->Draw("p0:iteration","fp0==0 && layer==10","same*l"); | |
558 | col(nteva,2); nteva->Draw("p0:iteration","fp0==0 && layer==11","same*l"); | |
559 | col(nteva,3); nteva->Draw("p0:iteration","fp0==0 && layer==12","same*l"); | |
560 | col(nteva,4); nteva->Draw("p0:iteration","fp0==0 && layer==13","same*l"); | |
561 | c->cd(3); du1->Draw(); | |
562 | col(nteva,1); nteva->Draw("fp1:iteration","p0==0 && layer==10","same*l"); | |
563 | col(nteva,2); nteva->Draw("fp1:iteration","p0==0 && layer==11","same*l"); | |
564 | col(nteva,3); nteva->Draw("fp1:iteration","p0==0 && layer==12","same*l"); | |
565 | col(nteva,4); nteva->Draw("fp1:iteration","p0==0 && layer==13","same*l"); | |
566 | c->cd(4); du2->Draw(); | |
567 | col(nteva,1); nteva->Draw("fp2:iteration","p0==0 && layer==10","same*l"); | |
568 | col(nteva,2); nteva->Draw("fp2:iteration","p0==0 && layer==11","same*l"); | |
569 | col(nteva,3); nteva->Draw("fp2:iteration","p0==0 && layer==12","same*l"); | |
570 | col(nteva,4); nteva->Draw("fp2:iteration","p0==0 && layer==13","same*l"); | |
571 | c->Print("amlsp.gif"); | |
572 | } | |
573 | //============================================================================= | |
574 | void align_module_layers_print(int iterations) { | |
575 | ||
576 | // print the result of align_module_layers | |
577 | ||
578 | int flayer = AliGeomManager::kTRD1; | |
579 | int llayer = AliGeomManager::kTRD6; | |
580 | for (int i=0; i<iterations; i++) { | |
581 | printf("\niteration %d\n",i); | |
582 | for (int layer=flayer; layer<=llayer; layer++) { | |
583 | printf("layer %2d ",layer); | |
584 | for (int j=0; j<6; j++) printf("%10.3f ",aml[i][layer][j]); | |
585 | printf("\n"); | |
586 | } | |
587 | } | |
588 | } | |
589 | //============================================================================= | |
590 | void align_totpc(char *inpfil, int what, int fitter_flag, int minim_flag, char *outfil) { | |
591 | ||
592 | // loop through detectors and align | |
593 | // track points will be read from inpfil??.root | |
594 | // what=0 ... detectors | |
595 | // what=1 ... stacks | |
596 | // what=2 ... supermodules | |
597 | // result will be stored on outfil.dat and outfil.stat | |
598 | ||
599 | AliTRDalignment ch; // misalignment in term of chambers | |
600 | AliTRDalignment sm; // misalignment in terms of supermodule (applies if what==3) | |
601 | int smpat[18] = {1,1,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,1}; | |
602 | int ntracks[540] = {}; | |
603 | enum {kDetector,kStack,kSupermodule}; | |
604 | for (int sec=0; sec<18; sec++) { | |
605 | if (!smpat[sec]) continue; | |
606 | TString filnam = Form("%s%02d.root",inpfil,sec); | |
607 | init_alt(filnam.Data(), fitter_flag, minim_flag); | |
608 | double pa[6]; | |
609 | int ntra=0; | |
610 | if (what==kSupermodule) { // align supermodels | |
611 | ntra = align_sm(sec, 1,0.1,pa,1); | |
612 | // use the mean between the two central chambers for the whole supermodule | |
613 | double pa0[6]; | |
614 | double pa1[6]; | |
615 | alt->GetAlignObj(AliGeomManager::LayerToVolUID(11,sec*5+2))->GetLocalPars(pa0,pa0+3); | |
616 | alt->GetAlignObj(AliGeomManager::LayerToVolUID(12,sec*5+2))->GetLocalPars(pa1,pa1+3); | |
617 | for (int i=0; i<6; i++) pa[i] = (pa0[i]+pa1[i])/2.0; | |
618 | sm.SetSm(sec,pa); | |
619 | } | |
620 | for (int st=0; st<5; st++) { | |
621 | int module = sec*5+st; | |
622 | if (what==kStack) ntra = align_volume(module, 9,14, 7,8, 1,0.1,pa,2); // align stacks | |
623 | for (int lay=0; lay<6; lay++) { | |
624 | if (sec==8 && st==4 && lay==3) continue; | |
625 | if (what==kDetector) ntra = align_volume(module, 9+lay,9+lay, 7,8, 1,0.1,pa,0); | |
626 | int det = AliTRDgeometry::GetDetector(lay,st,sec); | |
627 | ntracks[det] = ntra; | |
628 | alt->GetAlignObj(AliGeomManager::LayerToVolUID(9+lay,module))->GetLocalPars(pa,pa+3); // alternative way | |
629 | ch.SetCh(det,pa); | |
630 | // getchar(); | |
631 | } | |
632 | } | |
633 | } | |
634 | if (what==kSupermodule) sm.WriteAscii(Form("%ssm.dat",outfil)); | |
635 | ch.WriteAscii(Form("%s.dat",outfil)); | |
636 | FILE *fp=fopen(Form("%s.stat",outfil),"w"); | |
637 | for (int i=0; i<540; i++) fprintf(fp,"%3d %12d\n",i,ntracks[i]); | |
638 | fclose(fp); | |
639 | } | |
640 | //============================================================================= |