next50 trigger mask in AliHLTGlobalEsdConverterComponent
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDalignVolume.C
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 //=============================================================================