Update of calibration task and replace AliTRDrunCalib.C by AddTaskTRDCalib.C (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDalignVolume.C
CommitLineData
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
31AliAlignmentTracks *alt;
32AliTrackFitter *trf;
33AliTrackResiduals *res;
34TNtuple *ntresi; // residuals after minimization
35TNtuple *ntalpa; // alignment params
36TNtuple *nteva; // results of minimization
37double aml[100][40][6]; // results from align_module_layers [iter][layer][phi,z,r,rphi,rz,rr]
38enum {kStraight, kRieman, kKalman}; // track fitters
39enum {kLinear, kFast, kChi2}; // minimizer
40
41//=============================================================================
42TNtuple *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//=============================================================================
80void 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//=============================================================================
159void 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//=============================================================================
166TCanvas *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//=============================================================================
191void 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//=============================================================================
235void 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//=============================================================================
248int 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//=============================================================================
366int 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//=============================================================================
398int 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//=============================================================================
418void 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//=============================================================================
436void 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//=============================================================================
498void 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//=============================================================================
535void col(TNtuple *nt, int i) {
536 nt->SetLineColor(i);
537 nt->SetMarkerColor(i);
538}
539//=============================================================================
540void 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//=============================================================================
574void 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//=============================================================================
590void 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//=============================================================================