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 //=============================================================================