1 #if !defined(__CINT__) || defined(__MAKECINT__)
5 #include "TGeoManager.h"
15 #include "AliTRDgeometry.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"
29 //#include "/misc/misko/root/peakfit.C"
31 AliAlignmentTracks *alt;
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
41 //=============================================================================
42 TNtuple *makeNtuple(AliTrackResiduals *res, char *name, int flag) {
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
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
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();
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();
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);
79 //=============================================================================
80 void filterTracks(AliTrackResiduals *res, double level) {
82 // find the peak in y-residuals; remove outlier tracks
84 // first, find the peak and prepare the track quality cut
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");
92 if (hist->GetEntries()<100) return;
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;
100 // find the edges of the distribution
102 {for (bleft=bpeak; bleft>0; bleft--) {
103 if (hist->GetBinContent(bleft)<level*ypeak) break;
105 {for (brigt=bpeak; brigt<nbins; brigt++) {
106 if (hist->GetBinContent(brigt)<level*ypeak) break;
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);
116 // second, make copy of original tracks
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);
130 // third, fire all employees and hire only the good ones
132 res->SetNTracks(ntracks);
133 {for (Int_t itrack = 0; itrack < ntracks; itrack++) {
134 AliTrackPointArray *volarray = tmpvol[itrack];
135 AliTrackPointArray *traarray = tmptra[itrack];
137 {for (Int_t ipoint = 0; ipoint < volarray->GetNPoints(); ipoint++) {
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();
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;
155 if (good_track) res->AddTrackPointArrays(volarray, traarray);
158 //=============================================================================
159 void drawWithTitles(TTree *tr, char *what, char *xtitle, char *ytitle) {
161 tr->GetHistogram()->GetXaxis()->SetTitle(xtitle);
162 tr->GetHistogram()->GetYaxis()->SetTitle(ytitle);
163 tr->GetHistogram()->DrawCopy();
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();
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
185 // peakfit(nt->GetHistogram(),0.1,0.1,par);
186 nt->GetHistogram()->Fit("gaus");
187 ((TH1F*) gDirectory->Get("dyhist"))->DrawCopy();
190 //=============================================================================
191 void init_alt(const char *trapo, int fitter_flag, int res_flag) {
193 // prepare the AliAlignmentTracks object
194 // read space points from file trapo
196 alt = new AliAlignmentTracks();
197 alt->SetPointsFilename(trapo);
198 //alt->Misalign("kuku.root","TRDAlignObjs");
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);
209 trf->SetMinNPoints(30); // if working with clusters
210 trf->SetMinNPoints(4); // if working with tracklets
211 alt->SetTrackFitter(trf);
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
227 res->SetMinNPoints(10); // if working with clusters
228 res->SetMinNPoints(1); // if working with tracklets
229 alt->SetMinimizer(res);
231 if (nteva) delete nteva;
232 nteva = new TNtuple("eva","eva","chamber:module:layer:iteration:p0:p1:p2:p3:p4:p5:fp0:fp1:fp2");
234 //=============================================================================
235 void misa_alt(UShort_t volid, Double_t *mis) {
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.
244 AliAlignObjParams *al = (AliAlignObjParams *) alt->GetAlignObj(volid);
245 al->SetLocalPars(mis[0],mis[1],mis[2],mis[3],mis[4],mis[5]);
247 //=============================================================================
248 int align_volume(TArrayI volIds, TArrayI volIdRefs, int iterations,
249 double level, double *pa, int showflag) {
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
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)));
272 int resu = alt->AlignVolumes(&volIds, &volIdRefs, AliGeomManager::kTPC1, AliGeomManager::kTOF,iterations);
275 // repeat taking only peak
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
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));
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());
300 {for (int i=0; i<volIds.GetSize(); i++) {
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]);
312 if (ntresi) delete ntresi;
313 ntresi = makeNtuple(res,"kuku",showflag-1);
314 TCanvas *c = show(ntresi);
320 te.SetTextSizePixels(11);
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;
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);
342 te.DrawTextNDC(xpos,ypos-=0.07,"Result");
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");
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);
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");
363 return res->GetNFilledTracks();
365 //=============================================================================
366 int align_volume(int module, int layer0, int layer1, int reflayer0, int reflayer1,
367 int iterations, double level, double *pa, int showflag) {
369 // align the chambers (module,layer0-layer1) to the layers reflayer0-reflayer1
370 // (but excluding the layers to be aligned)
372 // prepare volume arrays
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);
384 if (i>=AliGeomManager::kTRD1 && i<=AliGeomManager::kTRD6) {
385 temp[n++] = AliGeomManager::LayerToVolUID(i,module);
388 TArrayI volIdRefs(n);
389 for (int i=0; i<n; i++) volIdRefs.AddAt(temp[i],i);
391 TArrayI volIds(layer1-layer0+1);
393 for (int i=layer0; i<=layer1; i++) volIds.AddAt(AliGeomManager::LayerToVolUID(i,module),i-layer0);
395 return align_volume(volIds,volIdRefs,iterations,level,pa,showflag);
397 //=============================================================================
398 int align_sm(int sec, int iterations, double level, double *pa, int showflag) {
400 // align the whole supermodule to the TPC
401 // could be unified with the align_volume above
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);
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++);
415 return align_volume(volIds,volIdRefs,iterations,level,pa,showflag);
417 //=============================================================================
418 void align_layer_modules(int layer, int reflayer0, int reflayer1,
419 int iterations, double level) {
421 // Loop through modules (TRD chambers) within given layer.
422 // Align each module using other layers.
423 // Store result on an ntuple.
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]);
435 //=============================================================================
436 void align_module_layers(int module, int iterations, double level) {
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.
443 int flayer = AliGeomManager::kTRD1;
444 int llayer = AliGeomManager::kTRD6;
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],
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);
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]);
472 if (ntresi) delete ntresi;
473 ntresi = makeNtuple(res,"kuku",1);
474 sprintf(buf,"dy>>hidy%02d_%d_afte(100,-2,2)",i,layer);
476 sprintf(buf,"hidy%02d_%d_afte",i,layer);
477 dy=(TH1F*) gDirectory->Get(buf);
478 // peakfit(dy,0.1,0.1,par);
480 nteva->Fill(det,module,layer,i+0.2,0,0,0,0,0,0,
481 par[0],par[1],par[2]);
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(
497 //=============================================================================
498 void align_module_layers_plot(int iterations, double ymax) {
499 // show histograms produced by align_module_layers
501 int flayer = AliGeomManager::kTRD1;
502 int llayer = AliGeomManager::kTRD6;
504 if (2*iterations>ncol) ncol=2*iterations;
505 TCanvas *c = new TCanvas();
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, "");
514 for (int layer=flayer; layer<=llayer; layer++) {
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);
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++;
534 //=============================================================================
535 void col(TNtuple *nt, int i) {
537 nt->SetMarkerColor(i);
539 //=============================================================================
540 void align_module_layers_summary_plot() {
544 TCanvas *c = new TCanvas("c","c");
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)");
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");
573 //=============================================================================
574 void align_module_layers_print(int iterations) {
576 // print the result of align_module_layers
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]);
589 //=============================================================================
590 void align_totpc(char *inpfil, int what, int fitter_flag, int minim_flag, char *outfil) {
592 // loop through detectors and align
593 // track points will be read from inpfil??.root
594 // what=0 ... detectors
596 // what=2 ... supermodules
597 // result will be stored on outfil.dat and outfil.stat
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);
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
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;
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);
628 alt->GetAlignObj(AliGeomManager::LayerToVolUID(9+lay,module))->GetLocalPars(pa,pa+3); // alternative way
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]);
640 //=============================================================================