ITS/AliITSRealignTracks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 #include <TArray.h>
  18 #include <TFile.h>
  19 #include <TStopwatch.h>
  20 #include <TArray.h>
  21 #include <TNtuple.h>
  22 #include <TClonesArray.h>
  23 #include <TMath.h>
  24 #include <TGeoManager.h>
  25 #include <TSystem.h>
  26 #include <TGeoMatrix.h>
  27 #include <TGraph.h>
  28 #include <TCanvas.h>
  29 #include <TH1F.h>
  30 #include "AliITSRealignTracks.h"
  31 #include "AliAlignmentTracks.h"
  32 #include "AliAlignObjParams.h"
  33 #include "AliAlignObj.h"
  34 #include "AliGeomManager.h"
  35 #include "AliTrackFitter.h"
  36 #include "AliTrackFitterKalman.h"
  37 #include "AliTrackFitterRieman.h"
  38 #include "AliTrackResidualsFast.h"
  39 #include "AliTrackResidualsChi2.h"
  40 #include "AliTrackResidualsLinear.h"
  41 #include "AliLog.h"
  42
  43 /* $Id$ */
  44
  45
  46 ClassImp(AliITSRealignTracks)
  47
  48 const Int_t referSect=2;
  49
  50
  51 AliITSRealignTracks::AliITSRealignTracks(TString minimizer,Int_t fit,Bool_t covUsed,TString fileintro,TString geometryfile,TString misalignmentFile,TString startingfile):
  52   AliAlignmentTracks(),
  53   fSurveyObjs(0),
  54   fgeomfilename(),
  55   fmintracks(),
  56   fCovIsUsed(covUsed),
  57   fUpdateCov(kFALSE),
  58   fVarySigmaY(kFALSE),
  59   fCorrModules(0),
  60   fLimitCorr(0),
  61   fsigmaY(),
  62   fDraw(kFALSE),
  63   fAlignDrawObjs(0),
  64   fCanvPar(0),
  65   fCanvGr(0),
  66   fgrIterMeanX(0),
  67   fgrIterRMSX(0),
  68   fgrIterMeanY(0),
  69   fgrIterRMSY(0),
  70   fgrIterMeanZ(0),
  71   fgrIterRMSZ(0),
  72   fgrIterMeanPsi(0),
  73   fgrIterRMSPsi(0),
  74   fgrIterMeanTheta(0),
  75   fgrIterRMSTheta(0),
  76   fgrIterMeanPhi(0),
  77   fgrIterRMSPhi(0)
  78 {
  79
  80   // minimizer="fast"->AliTrackResidualFast minimizer
  81   //           "minuit"->AliTrackResidualChi2 minimizer
  82   //           "minuitnorot"->AliTrackResidualChi2 minimizer without rotations degrees of freedom
  83   //           "linear"->AliTrackResidualLinear minimizer
  84   //fit=0-> Riemann Fitter, fit=1->Kalman
  85   //fileintro=file into which the Tree with the space points is stored
  86   //geometryfile=file containing the geometry
  87
  88
  89   SetPointsFilename(fileintro.Data());
  90   SetGeomFilename(geometryfile);
  91   InitAlignObjs();
  92   if(!InitSurveyObjs(kFALSE))AliWarning("Unable to set Survey AlignObjs!");
  93
  94   if(startingfile=="")printf("Starting from default geometry \n");
  95   else ReadAlignObjs(startingfile.Data(),"ITSAlignObjs");
  96
  97   if(misalignmentFile=="")printf("NO FAKE MISALIGNMENT INTRODUCED \n");
  98   else {
  99     Bool_t misal=Misalign(misalignmentFile,"ITSAlignObjs");
 100     if(!misal)AliWarning("Incorrect fake misalignment filename!");;
 101   }
 102
 103   if(!gGeoManager) AliGeomManager::LoadGeometry(fgeomfilename.Data());
 104   if(covUsed)SetCovIsUsed(kTRUE);
 105   if(!SelectFitter(fit))AliWarning("Incorrect fitter assignment!");
 106   if(!SelectMinimizer(minimizer))AliWarning("Incorrect minimizer assignment!");
 107   fsigmaY=1.;
 108   fmintracks=1;
 109   BuildIndex();
 110
 111 }
 112
 113 AliITSRealignTracks::AliITSRealignTracks(const AliITSRealignTracks &realignTracks):
 114   AliAlignmentTracks(),
 115   fSurveyObjs(new AliAlignObj**(*realignTracks.fSurveyObjs)),
 116   fgeomfilename(realignTracks.fgeomfilename),
 117   fmintracks(realignTracks.fmintracks),
 118   fCovIsUsed(realignTracks.fCovIsUsed),
 119   fUpdateCov(realignTracks.fUpdateCov),
 120   fVarySigmaY(realignTracks.fVarySigmaY),
 121   fCorrModules(new Double_t *(*realignTracks.fCorrModules)),
 122   fLimitCorr(realignTracks.fLimitCorr),
 123   fsigmaY(realignTracks.fsigmaY),
 124   fDraw(kFALSE),
 125   fAlignDrawObjs(realignTracks.fAlignDrawObjs),
 126   fCanvPar(realignTracks.fCanvPar),
 127   fCanvGr(realignTracks.fCanvGr),
 128   fgrIterMeanX(realignTracks.fgrIterMeanX),
 129   fgrIterRMSX(realignTracks.fgrIterRMSX),
 130   fgrIterMeanY(realignTracks.fgrIterMeanY),
 131   fgrIterRMSY(realignTracks.fgrIterRMSY),
 132   fgrIterMeanZ(realignTracks.fgrIterMeanZ),
 133   fgrIterRMSZ(realignTracks.fgrIterRMSZ),
 134   fgrIterMeanPsi(realignTracks.fgrIterMeanPsi),
 135   fgrIterRMSPsi(realignTracks.fgrIterRMSPsi),
 136   fgrIterMeanTheta(realignTracks.fgrIterMeanTheta),
 137   fgrIterRMSTheta(realignTracks.fgrIterRMSTheta),
 138   fgrIterMeanPhi(realignTracks.fgrIterMeanPhi),
 139   fgrIterRMSPhi(realignTracks.fgrIterRMSPhi)
 140
 141 {//Copy Constructor
 142   AliWarning("Can't copy AliAlignmentTracks Data member!");
 143 }
 144
 145 AliITSRealignTracks& AliITSRealignTracks::operator=(const AliITSRealignTracks &obj){
 146   ////////////////////////
 147   // Assignment operator
 148   ////////////////////////
 149   this->~AliITSRealignTracks();
 150   new(this) AliITSRealignTracks(obj);
 151   return *this;
 152 }
 153
 154 AliITSRealignTracks::~AliITSRealignTracks(){
 155   //destructor
 156
 157   if(fSurveyObjs)   DeleteSurveyObjs();
 158   if(fAlignDrawObjs) DeleteDrawHists();
 159   //delete [] fSurveyObjs;
 160
 161 }
 162
 163
 164 //_____________________________
 165 Bool_t AliITSRealignTracks::SelectFitter(Int_t fit,Int_t minTrackPoint){
 166
 167   if(fit==1){
 168      AliTrackFitterKalman *fitter= new AliTrackFitterKalman();
 169      fitter->SetMinNPoints(minTrackPoint);
 170      SetTrackFitter(fitter);
 171   }
 172
 173   else if(fit==0){
 174     AliTrackFitterRieman *fitter=new AliTrackFitterRieman();
 175     fitter=new AliTrackFitterRieman();
 176     fitter->SetMinNPoints(minTrackPoint);
 177     SetTrackFitter(fitter);
 178   }
 179   else return kFALSE;
 180
 181   return kTRUE;
 182 }
 183
 184
 185 Bool_t AliITSRealignTracks::SelectMinimizer(TString minimizer,Int_t minpoints,const Bool_t *coord){
 186   AliTrackResiduals *res;
 187   if(minimizer=="minuit"){
 188     res = new AliTrackResidualsChi2();
 189     if(coord){
 190       for(Int_t j=0;j<6;j++){
 191         if(coord[j])res->FixParameter(j);
 192       }
 193     }
 194   }
 195   else if(minimizer=="minuitnorot"){
 196     res = new AliTrackResidualsChi2();
 197     res->FixParameter(3);
 198     res->FixParameter(4);
 199     res->FixParameter(5);
 200   }
 201
 202   else if(minimizer=="fast"){
 203     res = new AliTrackResidualsFast();
 204     if(coord){
 205       for(Int_t j=0;j<6;j++){
 206         if(coord[j])res->FixParameter(j);
 207       }
 208     }
 209   }
 210   else if(minimizer=="linear"){
 211     res = new AliTrackResidualsLinear();
 212   }
 213
 214   else {
 215     printf("Trying to set a non existing minimizer! \n");
 216     return kFALSE;
 217   }
 218
 219   res->SetMinNPoints(minpoints);
 220   SetMinimizer(res);
 221
 222   return kTRUE;
 223 }
 224
 225 //____________________________________
 226 void AliITSRealignTracks::SetVarySigmaY(Bool_t varysigmay,Double_t sigmaYfixed){
 227
 228   fVarySigmaY=varysigmay;
 229   if(!varysigmay){
 230     if(sigmaYfixed>0.)fsigmaY=sigmaYfixed;
 231     else {
 232       printf("Negative assignment to sigmaY! set it to default value (1 cm) \n");
 233       fsigmaY=1.;
 234     }
 235   }
 236 }
 237
 238 //_______________________________________
 239 void AliITSRealignTracks::RealignITSVolIndependent(Int_t iter1,Int_t iterations,Int_t minNtracks,Int_t layer,Int_t minTrackPoint){
 240
 241
 242   //iter1=#iterations inside AliAlignmentTracks::AliAlignVolumesITS method
 243   //iterations=#iterations on the all procedure
 244   //layer=0->all ITS, otherways the usual notation is considered (1=SPD1,2=SPD2,3=SDD1,4=SDD2,5=SSD1,6=SSD2)
 245   //minNtracks=minimun number of tracks passing through a module in order to try to realign the module itsself
 246   //           if minNtracks<0, minimun number of tracks is |minNtracks|*minNumPoint[layer]/fact (see the code below): this allows a different
 247   //           choice of the number of tracks required on different layers and to vary these numbers once tuned the relative proportions.
 248   //minTrackPoint=minimun number of "good" points required to a track (THE POINT ON THE MODULE THAT IS GOING TO BE REALIGNED
 249   //IS NEVER CONSIDERED->max number can be required is 11 for cosmics tracks) for the track being considered in the minimization
 250
 251
 252   fTrackFitter->SetMinNPoints(minTrackPoint);
 253   TArrayI volIDs2(2200);
 254   volIDs2.Reset(0);
 255   TArrayI volIDs(1);
 256   TString command;
 257   TArrayI volIDsFit;
 258
 259   Int_t iLayer,iLayerToAlign;
 260
 261   Int_t minNumPoint[6]={100,100,100,100,50,50};
 262   Double_t fact=10;
 263   Int_t j=0;
 264
 265   Int_t size=0;
 266   Int_t layerNum,modNum,lastVolid=0;
 267   TNtuple *ntVolumeAlign=new TNtuple("ntVolumeAlign","NTuple with volume tried to be realigned","layerNum:modNum:volumeIDnum");
 268
 269   TStopwatch *timer=new TStopwatch();
 270   timer->Start();
 271   BuildIndex();
 272   j=0;
 273   UShort_t volid;
 274
 275   for(Int_t iter=0;iter<iterations;iter++){
 276
 277     //Starting Independent Modules Realignment
 278     for(iLayerToAlign=(Int_t)AliGeomManager::kSPD1;iLayerToAlign<=(Int_t)AliGeomManager::kSSD2;iLayerToAlign++){
 279       if(layer!=0&&iLayerToAlign!=layer)continue;
 280       j=0;
 281       size=0;
 282       for(Int_t k=(Int_t)AliGeomManager::kSPD1;k<=(Int_t)AliGeomManager::kSSD2;k++){
 283         size+=AliGeomManager::LayerSize(k);
 284         printf("size: %d \n",size);
 285       }
 286
 287       for (Int_t iModule=0;iModule<AliGeomManager::LayerSize(iLayerToAlign);iModule++){
 288         j=0;
 289         if(minNtracks<0){
 290           if(GetLastIndex(iLayerToAlign-(Int_t)AliGeomManager::kFirstLayer,iModule)<minNumPoint[iLayerToAlign-(Int_t)AliGeomManager::kFirstLayer]*(-1*minNtracks/fact))continue;        }
 291         else if(GetLastIndex(iLayerToAlign-(Int_t)AliGeomManager::kFirstLayer,iModule)<minNtracks)continue;
 292
 293         UShort_t volidAl = AliGeomManager::LayerToVolUID(iLayerToAlign,iModule);
 294
 295         volIDsFit.Reset(0);
 296         volIDsFit.Set(size-1);
 297         for (iLayer=AliGeomManager::kSPD1;iLayer<AliGeomManager::kTPC1;iLayer++){
 298           for (Int_t iModule2=0;iModule2<AliGeomManager::LayerSize(iLayer);iModule2++){
 299             volid = AliGeomManager::LayerToVolUID(iLayer,iModule2);
 300             if(AliGeomManager::LayerToVolUID(iLayer,iModule2)==volidAl)continue;
 301             volIDsFit.AddAt(volid,j);
 302             j++;
 303           }
 304         }
 305         volIDs.AddAt((Int_t)volidAl,0);
 306         if(iter==iterations-1){
 307           volIDs2.AddAt(volidAl,lastVolid);
 308           lastVolid++;
 309         }
 310         volIDs2.AddAt(volidAl,lastVolid);
 311         AlignVolumesITS(&volIDs,&volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,iter1);
 312       }
 313     }
 314
 315
 316     if((iter+1)%5==0||iter==0||iter==1||iter==2||iter==3||iter==iterations-1){
 317       command="RealignObj";
 318       command+=(iter+1);
 319       command.Append(".root");
 320       WriteRealignObjArray(command.Data(),AliGeomManager::kSPD1,AliGeomManager::kSSD2);
 321     }
 322   }
 323
 324
 325   if(j==0){printf("j=0 \n");return;}
 326   for(Int_t k=0;k<volIDs2.GetSize();k++){
 327     if(volIDs2.At(k)==0)break;
 328     layerNum=AliGeomManager::VolUIDToLayer(volIDs2.At(k),modNum);
 329     ntVolumeAlign->Fill(layerNum,modNum,volIDs2.At(k));
 330   }
 331   printf("End of selecting modules cycle: %d modules selected \n",j);
 332   TFile *f=new TFile("RealignVolNt.root","RECREATE");
 333   f->cd();
 334   ntVolumeAlign->Write();
 335   f->Close();
 336
 337   timer->Stop();
 338   timer->Print();
 339   return;
 340 }
 341
 342
 343 void AliITSRealignTracks::RealignITStracks(TString minimizer,Int_t fit=0,Int_t iter1=1,Int_t iterations=5,Int_t minNtracks=-10,Int_t layer=0,Int_t minTrackPoint=6,Bool_t covUsed=kFALSE,TString misalignmentFile="",TString startingfile="",Int_t doGlobal=1)
 344 {
 345
 346   // minimizer="fast"->AliTrackResidualFast minimizer
 347   //           "minuit"->AliTrackResidualChi2 minimizer
 348   //           "minuitnorot"->AliTrackResidualChi2 minimizer without rotations degrees of freedom
 349   //           "linear"->AliTrackResidualLinear minimizer
 350   //fit=0-> Riemann Fitter, fit=1->Kalman
 351   //iter1=#iterations inside AliAlignmentTracks::AliAlignVolumesITS method
 352   //iterations=#iterations on the all procedure
 353   //layer=0->all ITS, otherways the usual notation is considered (1=SPD1,2=SPD2,3=SDD1,4=SDD2,5=SSD1,6=SSD2)
 354   //minNtracks=minimun number of tracks passing through a module in order to try to realign the module itsself
 355   //           if minNtracks<0, minimun number of tracks is |minNtracks|*minNumPoint[layer]/fact (see the code below): this allows a different
 356   //           choice of the number of tracks required on different layers and to vary these numbers once tuned the relative proportions.
 357   //minTrackPoint=minimun number of "good" points required to a track (THE POINT ON THE MODULE THAT IS GOING TO BE REALIGNED
 358   //              IS NEVER CONSIDERED->max number that can be required is 11 for cosmics tracks) for the track being considered in the minimization
 359   //doGlobal : do global realignment, 0=no, 1= yes, 2=only global
 360
 361
 362   TArrayI volIDs2(2200);
 363   volIDs2.Reset(0);
 364   TArrayI volIDs(1);
 365   TString command;
 366   TArrayI volIDsFit;
 367
 368   Int_t iLayer,iLayerToAlign;
 369
 370   Int_t minNumPoint[6]={100,100,100,100,50,50};
 371   Double_t fact=10;
 372   Int_t count=0;
 373
 374   Int_t size=0;
 375   Int_t layerNum,modNum,lastVolid=0;
 376   TNtuple *ntVolumeAlign=new TNtuple("ntVolumeAlign","NTuple with volume tried to be realigned","layerNum:modNum:volumeIDnum");
 377
 378
 379   if(!SelectFitter(fit))AliWarning("Incorrect fitter assignment!");
 380   if(!SelectMinimizer(minimizer))AliWarning("Incorrect minimizer assignment!");
 381   if(misalignmentFile=="")printf("NO FAKE MISALIGNMENT INTRODUCED \n");
 382   else {
 383     Bool_t misal=Misalign(misalignmentFile,"ITSAlignObjs");
 384     if(!misal)return;
 385   }
 386
 387
 388   TStopwatch *timer=new TStopwatch();
 389   timer->Start();
 390   BuildIndex();
 391   count=0;
 392   UShort_t volid;
 393
 394   if(startingfile=="")printf("Starting from default geometry \n");
 395   else {
 396     printf("Starting from AlignObjs file: %s",startingfile.Data());
 397     ReadAlignObjs(startingfile.Data(),"ITSAlignObjs");
 398   }
 399
 400   for(Int_t iter=0;iter<iterations;iter++){
 401     if(covUsed)SetCovIsUsed(kTRUE);
 402
 403
 404     //START HIERARCHY REALIGNMENT
 405
 406     if(layer==0&&(doGlobal==1||doGlobal==2)){
 407       for(Int_t siter=0;siter<5;siter++){
 408         fTrackFitter->SetMinNPoints(2);
 409         SetCovUpdate(kFALSE);
 410         AlignSPDHalfBarrelToSectorRef(referSect,3);
 411         //      AlignSPDBarrel(1);
 412         //      if(siter==0)SetCovUpdate(kFALSE);
 413         //      AlignSPDHalfBarrel(0,3);
 414         //      SetCovUpdate(kTRUE);
 415         AlignSPDHalfBarrelToHalfBarrel(1,3);
 416         //      AlignSPDHalfBarrelToSectorRef(referSect,3);
 417         for(Int_t sector=0;sector<10;sector++){
 418           SetMinNtracks(100);
 419           if(sector==referSect)continue;
 420           AlignSPDSectorWithSectors(sector,1);
 421         }
 422
 423
 424         for(Int_t lay=1;lay<=6;lay++){
 425           if(!AlignLayerToSector(lay,referSect,3))AlignLayerToSPDHalfBarrel(lay,0,3);
 426         }
 427         AlignSPDHalfBarrel(0,3);
 428
 429         Int_t layers[6]={2,2,1,0,0,0};
 430         fTrackFitter->SetMinNPoints(4);
 431         AlignLayersToLayers(layers,1);
 432
 433         fTrackFitter->SetMinNPoints(6);
 434         layers[2]=2;
 435         layers[3]=1;//{2,2,2,1,0,0};
 436         AlignLayersToLayers(layers,1);
 437
 438         fTrackFitter->SetMinNPoints(6);
 439         layers[3]=2;
 440         layers[4]=1;//{2,2,2,2,1,0};
 441         AlignLayersToLayers(layers,1);
 442
 443         fTrackFitter->SetMinNPoints(6);
 444         layers[4]=2;
 445         layers[5]=1;//{2,2,2,2,2,1};
 446         AlignLayersToLayers(layers,1);
 447
 448
 449         for(Int_t sector=0;sector<10;sector++){
 450           AlignSPDSectorToOuterLayers(sector,1);
 451         }
 452         WriteRealignObjArray("AfterGlobal.root",AliGeomManager::kSPD1,AliGeomManager::kSSD2);
 453       }
 454     }
 455
 456     if(doGlobal==2)return;
 457
 458     if(covUsed)SetCovUpdate(kTRUE);
 459     SetMinNtracks(1);
 460
 461
 462     // STARTS INDEPENDENT MOULES REALIGNMENT
 463
 464     fTrackFitter->SetMinNPoints(minTrackPoint);
 465     for(iLayerToAlign=(Int_t)AliGeomManager::kSPD1;iLayerToAlign<=(Int_t)AliGeomManager::kSSD2;iLayerToAlign++){
 466       if(layer!=0&&iLayerToAlign!=layer)continue;
 467       count=0;
 468       size=0;
 469       for(Int_t k=(Int_t)AliGeomManager::kSPD1;k<=(Int_t)AliGeomManager::kSSD2;k++){
 470         size+=AliGeomManager::LayerSize(k);
 471         printf("size: %d \n",size);
 472       }
 473
 474       for (Int_t iModule=0;iModule<AliGeomManager::LayerSize(iLayerToAlign);iModule++){
 475         count=0;
 476         if(minNtracks<0){
 477           if(GetLastIndex(iLayerToAlign-(Int_t)AliGeomManager::kFirstLayer,iModule)<minNumPoint[iLayerToAlign-(Int_t)AliGeomManager::kFirstLayer]*(-1*minNtracks/fact))continue;
 478         }
 479         else if(GetLastIndex(iLayerToAlign-(Int_t)AliGeomManager::kFirstLayer,iModule)<minNtracks)continue;
 480
 481         UShort_t volidAl = AliGeomManager::LayerToVolUID(iLayerToAlign,iModule);
 482
 483         volIDsFit.Reset(0);
 484         volIDsFit.Set(size-1);
 485         for (iLayer=(Int_t)AliGeomManager::kSPD1;iLayer<(Int_t)AliGeomManager::kTPC1;iLayer++){
 486           for (Int_t iModule2=0;iModule2<AliGeomManager::LayerSize(iLayer);iModule2++){
 487
 488
 489             volid = AliGeomManager::LayerToVolUID(iLayer,iModule2);
 490
 491             if(AliGeomManager::LayerToVolUID(iLayer,iModule2)==volidAl)continue;
 492             volIDsFit.AddAt(volid,count);
 493             count++;
 494           }
 495         }
 496
 497         volIDs.AddAt((Int_t)volidAl,0);
 498         if(iter==iterations-1){
 499           volIDs2.AddAt(volidAl,lastVolid);
 500           lastVolid++;
 501         }
 502         volIDs2.AddAt(volidAl,lastVolid);
 503         AlignVolumesITS(&volIDs,&volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,iter1);
 504       }
 505     }
 506
 507     if((iter+1)%2==0||(iter+1)%5==0||iter==0||iter==1||iter==2||iter==3||iter==iterations-1){
 508       command="RealignObj";
 509       command+=(iter+1);
 510       command.Append(".root");
 511       WriteRealignObjArray(command.Data(),AliGeomManager::kSPD1,AliGeomManager::kSSD2);
 512     }
 513
 514   }
 515
 516   if(count==0){printf("count=0 \n");return;}
 517   for(Int_t k=0;k<volIDs2.GetSize();k++){
 518     if(volIDs2.At(k)==0)break;
 519     layerNum=AliGeomManager::VolUIDToLayer(volIDs2.At(k),modNum);
 520     ntVolumeAlign->Fill(layerNum,modNum,volIDs2.At(k));
 521   }
 522   printf("End of selecting modules cycle: %d modules selected \n",count);
 523   TFile *f=new TFile("RealignVolNt.root","RECREATE");
 524   f->cd();
 525   ntVolumeAlign->Write();
 526   f->Close();
 527
 528   timer->Stop();
 529   timer->Print();
 530   return;
 531
 532 }
 533
 534
 535 //______________________________________________________________________________
 536 void AliITSRealignTracks::InitAlignObjs()
 537 {
 538   // Initialize the alignment objects array
 539   TMatrixDSym c(6);
 540   Double_t cov[21];
 541   for(Int_t i=0;i<21;i++)cov[i]=0.;
 542   for(Int_t i=0;i<3;i++)cov[i*(i+1)/2+i]=0.05*0.05;//Set Default Error to 500 micron for Translations
 543   for(Int_t i=3;i<6;i++)cov[i*(i+1)/2+i]=0.001*0.001*180*180/3.14/3.14;//and 1 mrad for rotations (global ref. sysytem->~40 micron for SPD1,~450 micron for SSD2)
 544
 545   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 546   fAlignObjs = new AliAlignObj**[nLayers];
 547   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 548     fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 549     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 550       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
 551       fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
 552       fAlignObjs[iLayer][iModule]->SetCorrMatrix(cov);
 553       fAlignObjs[iLayer][iModule]->SetUniqueID(0);
 554     }
 555   }
 556 }
 557
 558 //______________________________________________________________________________
 559 void AliITSRealignTracks::ResetCorrModules(){
 560   // Initialize and reset to 0 the array with the information on correlation
 561   if(!fCorrModules){
 562     Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 563     fCorrModules = new Double_t*[nLayers];
 564     for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 565       fCorrModules[iLayer] = new Double_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 566       for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 567         fCorrModules[iLayer][iModule]=0.;
 568       }
 569     }
 570   }
 571   else{
 572     for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 573       for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 574         fCorrModules[iLayer][iModule]=0.;
 575       }
 576     }
 577   }
 578 }
 579
 580 //______________________________________________________________________________
 581 Bool_t AliITSRealignTracks::InitSurveyObjs(Bool_t infinite,Double_t factor,TString filename,TString arrayName){
 582
 583   if(fSurveyObjs)DeleteSurveyObjs();
 584   Bool_t fromfile=kFALSE;
 585   TFile *surveyObj;
 586   TClonesArray *clnarray;
 587   if(!filename.IsNull()){
 588     //Initialize from file
 589     if(gSystem->AccessPathName(filename.Data(),kFileExists)){
 590       printf("Wrong Survey AlignObjs File Name \n");
 591       return kFALSE;
 592     }
 593     if(arrayName.IsNull()){
 594       printf("Null Survey Object Name! \n");
 595       return kFALSE;
 596     }
 597
 598     fromfile=kTRUE;
 599   }
 600
 601   // Initialize the alignment objects array with default values
 602   Double_t v=1.*factor;
 603   if(infinite)v*=100000.;
 604   TMatrixDSym c(6);
 605   Double_t cov[21];
 606   for(Int_t i=0;i<21;i++)cov[i]=0.;
 607   for(Int_t i=0;i<3;i++)cov[i*(i+1)/2+i]=0.1*0.1*v;//Set Default Error to 1 mm  for Translation
 608   for(Int_t i=3;i<6;i++)cov[i*(i+1)/2+i]=0.03*0.03*180.*180./3.14/3.14*v;//and 30 mrad (~1.7 degrees)for rotations (global ref. sysytem)
 609
 610   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 611   fSurveyObjs = new AliAlignObj**[nLayers];
 612   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 613     fSurveyObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 614     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 615       fSurveyObjs[iLayer][iModule] = 0x0;
 616     }
 617   }
 618
 619   if(fromfile){
 620     surveyObj=TFile::Open(filename.Data());
 621     if (!surveyObj || !surveyObj->IsOpen()) {
 622       AliError(Form("Could not open SurveyObjs file: %s !",filename.Data()));
 623       return kFALSE;
 624     }
 625     printf("Getting TClonesArray \n");
 626     clnarray=(TClonesArray*)surveyObj->Get(arrayName);
 627     Int_t size=clnarray->GetSize();
 628     UShort_t volid;
 629     for(Int_t ivol=0;ivol<size;ivol++){
 630       AliAlignObjParams *a=(AliAlignObjParams*)clnarray->At(ivol);
 631       volid=a->GetVolUID();
 632       Int_t iModule;
 633       AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 634       if(iLayer<=0)continue;
 635       if(a->GetUniqueID()==0)continue;
 636       printf("Updating survey for volume: %d ,layer: %d module: %d from file\n",volid,iLayer,iModule);
 637       fSurveyObjs[iLayer-AliGeomManager::kFirstLayer][iModule] = new AliAlignObjParams(*a);
 638       fSurveyObjs[iLayer-AliGeomManager::kFirstLayer][iModule]->SetUniqueID(a->GetUniqueID());
 639       fSurveyObjs[iLayer-AliGeomManager::kFirstLayer][iModule]->Print("");
 640     }
 641     delete clnarray;
 642     surveyObj->Close();
 643   }
 644
 645   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 646     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 647       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
 648       if(!fSurveyObjs[iLayer][iModule]){
 649         printf("Updating survey for volume: %d ,layer: %d module: %d with default values \n",volid,iLayer,iModule);
 650         fSurveyObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
 651         fSurveyObjs[iLayer][iModule]->SetCorrMatrix(cov);
 652         fSurveyObjs[iLayer][iModule]->SetUniqueID(0);
 653       }
 654
 655     }
 656   }
 657
 658
 659   return kTRUE;
 660 }
 661
 662
 663 //______________________________________________________________________________
 664 Int_t AliITSRealignTracks::CheckWithSurvey(Double_t factor,TArrayI *volids){
 665
 666   // Check the parameters of the alignment objects in volids (or of all objects if volids is null)
 667   // are into the boundaries set by the cov. matrix of the survey objs
 668   // Returns the number of objects out of boudaries
 669   AliAlignObj *alignObj;
 670   Int_t outofsurv=0;
 671   UShort_t volid;
 672   Double_t surveycov[21],transl[3],rot[3],survtransl[3],survrot[3];
 673   if(volids==0x0){
 674     for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 675       for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++){
 676         volid=AliGeomManager::LayerToVolUIDSafe(iLayer+AliGeomManager::kFirstLayer,iModule);
 677         alignObj=GetAlignObj(volid);
 678         alignObj->GetPars(transl,rot);
 679         fSurveyObjs[iLayer][iModule]->GetCovMatrix(surveycov);
 680         fSurveyObjs[iLayer][iModule]->GetPars(survtransl,survrot);
 681         if(TMath::Sqrt(TMath::Abs(surveycov[0]))*factor<TMath::Abs(transl[0]-survtransl[0])||TMath::Sqrt(TMath::Abs(surveycov[2]))*factor<TMath::Abs(transl[1]-survtransl[1])||TMath::Sqrt(TMath::Abs(surveycov[5]))*factor<TMath::Abs(transl[2]-survtransl[2])||TMath::Sqrt(TMath::Abs(surveycov[9]))*factor<TMath::Abs(rot[0]-survrot[0])||TMath::Sqrt(TMath::Abs(surveycov[14]))*factor<TMath::Abs(rot[1]-survrot[1])||TMath::Sqrt(TMath::Abs(surveycov[20]))*factor<TMath::Abs(rot[2]-survrot[2])){
 682           printf("Results for module %d out of Survey: reinitializing it from survey \n",volid);
 683           //      *alignObj = *alignObjSurv;
 684           alignObj->SetPars(survtransl[0],survtransl[1],survtransl[2],survrot[0],survrot[1],survrot[2]);
 685           alignObj->SetUniqueID(0);
 686           if(fUpdateCov)alignObj->SetCorrMatrix(surveycov);
 687           outofsurv++;
 688         }
 689       }
 690     }
 691   }
 692   else{
 693     Int_t iLayer;
 694     Int_t iModule;
 695     for(Int_t j=0;j<volids->GetSize();j++){
 696       volid=volids->At(j);
 697       alignObj=GetAlignObj(volid);
 698       alignObj->GetPars(transl,rot);
 699       iLayer=(Int_t)AliGeomManager::VolUIDToLayerSafe(volid,iModule)-(Int_t)AliGeomManager::kFirstLayer;
 700       fSurveyObjs[iLayer][iModule]->GetCovMatrix(surveycov);
 701       fSurveyObjs[iLayer][iModule]->GetPars(survtransl,survrot);
 702       if(TMath::Sqrt(TMath::Abs(surveycov[0]))*factor<TMath::Abs(transl[0]-survtransl[0])||TMath::Sqrt(TMath::Abs(surveycov[2]))*factor<TMath::Abs(transl[1]-survtransl[1])||TMath::Sqrt(TMath::Abs(surveycov[5]))*factor<TMath::Abs(transl[2]-survtransl[2])||TMath::Sqrt(TMath::Abs(surveycov[9]))*factor<TMath::Abs(rot[0]-survrot[0])||TMath::Sqrt(TMath::Abs(surveycov[14]))*factor<TMath::Abs(rot[1]-survrot[1])||TMath::Sqrt(TMath::Abs(surveycov[20]))*factor<TMath::Abs(rot[2]-survrot[2])){
 703         printf("Results for module %d out of Survey: reinitializing it from survey \n",volid);
 704         //        *alignObj = *alignObjSurv;
 705         alignObj->SetPars(survtransl[0],survtransl[1],survtransl[2],survrot[0],survrot[1],survrot[2]);
 706         alignObj->SetUniqueID(0);
 707         if(fUpdateCov)alignObj->SetCorrMatrix(surveycov);
 708         outofsurv++;
 709       }
 710     }
 711   }
 712   return outofsurv;
 713 }
 714
 715 //___________________________________________________________________
 716
 717 void AliITSRealignTracks::ResetAlignObjs(Bool_t all,TArrayI *volids)
 718 {
 719   // Reset the alignment objects in volids or all if all=kTRUE
 720   if(all){
 721     for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 722       for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
 723         fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
 724     }
 725   }
 726   else{
 727     Int_t layer;
 728     Int_t mod;
 729     for(Int_t j=0;j<volids->GetSize();j++){
 730       layer=(Int_t)AliGeomManager::VolUIDToLayer(volids->At(j),mod)-(Int_t)AliGeomManager::kFirstLayer;
 731       fAlignObjs[layer][mod]->SetPars(0,0,0,0,0,0);
 732     }
 733   }
 734 }
 735
 736 //______________________________________________-
 737 void AliITSRealignTracks::DeleteSurveyObjs()
 738 {
 739   if(!fSurveyObjs)return;
 740   // Delete the alignment objects array
 741   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 742     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++){
 743       if (fSurveyObjs[iLayer][iModule]) delete fSurveyObjs[iLayer][iModule];
 744     }
 745
 746     if(fSurveyObjs[iLayer])delete [] fSurveyObjs[iLayer];
 747   }
 748
 749   delete [] fSurveyObjs;
 750   fSurveyObjs = 0;
 751 }
 752
 753
 754 //______________________________________________________________________________
 755 Bool_t AliITSRealignTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName){
 756
 757   // Read alignment object from a file: update the alignobj already present with the one in the file
 758   // To be replaced by a call to CDB
 759
 760   if(gSystem->AccessPathName(alignObjFileName,kFileExists)){
 761     printf("Wrong AlignObjs File Name \n");
 762     return kFALSE;
 763   }
 764
 765   TFile *fRealign=TFile::Open(alignObjFileName);
 766   if (!fRealign || !fRealign->IsOpen()) {
 767     AliError(Form("Could not open Align Obj File file %s !",alignObjFileName));
 768     return kFALSE;
 769   }
 770   printf("Getting TClonesArray \n");
 771   TClonesArray *clnarray=(TClonesArray*)fRealign->Get(arrayName);
 772   Int_t size=clnarray->GetSize();
 773   UShort_t volid;
 774
 775   for(Int_t ivol=0;ivol<size;ivol++){
 776     AliAlignObjParams *a=(AliAlignObjParams*)clnarray->At(ivol);
 777     volid=a->GetVolUID();
 778     Int_t iModule;
 779     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 780     if(iLayer<AliGeomManager::kFirstLayer||iLayer>AliGeomManager::kSSD2)continue;
 781     printf("Updating volume: %d ,layer: %d module: %d \n",volid,iLayer,iModule);
 782     *fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule] *= *a;
 783   }
 784
 785   delete clnarray;
 786   fRealign->Close();
 787   return kTRUE;
 788 }
 789
 790 //_________________________________________
 791 Bool_t AliITSRealignTracks::FirstAlignmentLayers(Bool_t *layers,Int_t minNtracks,Int_t iterations,Bool_t fitall,TArrayI *volidsSet){
 792
 793   //Align all modules in the set of layers independently according to a sequence based on the number of tracks passing through a given module
 794
 795   BuildIndex();
 796   TString name="DrawFirstAlignment_Layers";
 797   UShort_t voluid;
 798   Int_t **lastIndex;
 799   Int_t laymax = 0;
 800   Int_t modmax = 0;
 801   Int_t maxntr=0,nMod=0,modAligned=0,size=0;
 802   for(Int_t i=0;i<6;i++){
 803     if(layers[i]==1){
 804       size+=AliGeomManager::LayerSize(i+AliGeomManager::kFirstLayer);
 805       name+=i+1;
 806     }
 807   }
 808
 809   TArrayI *volFit=new TArrayI(size);
 810   TArrayI *volFit2=new TArrayI(size-1);
 811   TArrayI *sequence=new TArrayI(size);
 812   TArrayI *volIn=new TArrayI(1);
 813
 814   // Initialize the index arrays
 815   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 816   lastIndex = new Int_t*[nLayers];
 817   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 818     lastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 819     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 820       lastIndex[iLayer][iModule] =  fLastIndex[iLayer][iModule];
 821       if(iLayer<=(AliGeomManager::kSSD2-AliGeomManager::kFirstLayer)&&layers[iLayer]==1){
 822         volFit->AddAt(AliGeomManager::LayerToVolUID(iLayer+AliGeomManager::kFirstLayer,iModule),maxntr);
 823         maxntr++;
 824       }
 825     }
 826   }
 827   Int_t found=0;
 828   maxntr=minNtracks+1;
 829   while (maxntr>minNtracks){
 830     maxntr=minNtracks;
 831     for (Int_t iLayer = 0; iLayer <= (AliGeomManager::kSSD2 - AliGeomManager::kFirstLayer); iLayer++) {
 832       if(layers[iLayer]==0)continue;
 833       for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 834         if(lastIndex[iLayer][iModule]>maxntr){
 835           maxntr=lastIndex[iLayer][iModule];
 836           laymax=iLayer;
 837           modmax=iModule;
 838         }
 839       }
 840     }
 841     if(maxntr>minNtracks){
 842       voluid=AliGeomManager::LayerToVolUID(laymax+AliGeomManager::kFirstLayer,modmax);
 843       sequence->AddAt(voluid,nMod);
 844       lastIndex[laymax][modmax]=0;
 845       nMod++;
 846     }
 847   }
 848
 849   sequence->Set(nMod);
 850
 851  Int_t ilayer,imod;
 852   for(Int_t iter=0;iter<iterations;iter++){
 853     if(iter>0&&fDraw)UpdateDraw(sequence,iter,iter);
 854     modAligned=0;
 855     for(Int_t k=0;k<nMod;k++){
 856       TArrayI *volFit3;
 857       voluid=sequence->At(k);
 858       ilayer=AliGeomManager::VolUIDToLayer(voluid,imod);
 859       volIn->AddAt(voluid,0);
 860       found=0;
 861       if(!fitall){
 862         for(Int_t j=0;j<nMod;j++){
 863           if(j==k){
 864             found=1;
 865             continue;
 866           }
 867           else volFit2->AddAt(sequence->At(j),j-found);
 868         }
 869         volFit2->Set(nMod-1);
 870       }
 871       else{
 872         for(Int_t j=0;j<volFit->GetSize();j++){
 873           if(volFit->At(j)!=volIn->At(0))volFit2->AddAt(volFit->At(j),j-found);
 874           else found=1;
 875         }
 876       }
 877
 878       if(volidsSet){
 879         volFit3=IntersectVolArray(volidsSet,volFit2);
 880       }
 881       else volFit3=new TArrayI(*volFit2);
 882
 883
 884       if(AlignVolumesITS(volIn,volFit3,AliGeomManager::kSPD1,AliGeomManager::kTPC1,2))modAligned++;
 885       delete volFit3;
 886
 887     }
 888   }
 889   Int_t noutofsurv=CheckWithSurvey(2.,sequence);
 890   printf("%d modules into the sequence \n %d modules re-aligned \n %d modules moved far away from survey (-> reset) \n",nMod,modAligned,noutofsurv);
 891   name.Append("_iter");
 892   name+=iterations;
 893   name.Append(".root");
 894   if(fDraw)WriteHists(name.Data());
 895   delete volFit;
 896   delete volFit2;
 897   delete sequence;
 898
 899   return kTRUE;
 900
 901 }
 902
 903 //__________________________________________
 904 Bool_t AliITSRealignTracks::FirstAlignmentSPD(Int_t minNtracks,Int_t iterations,Bool_t fitall,TArrayI *volidsSet){
 905
 906   BuildIndex();
 907
 908   UShort_t voluid;
 909   Int_t **lastIndex;
 910   Int_t laymax = 0;
 911   Int_t modmax = 0;
 912   Int_t maxntr=0,nMod=0,modAligned=0;
 913   TArrayI *volFit=new TArrayI(AliGeomManager::LayerSize(1)+AliGeomManager::LayerSize(2));
 914   TArrayI *volFit2=new TArrayI(AliGeomManager::LayerSize(1)+AliGeomManager::LayerSize(2)-1);
 915   TArrayI *sequence=new TArrayI(AliGeomManager::LayerSize(1)+AliGeomManager::LayerSize(2));
 916   TArrayI *volIn=new TArrayI(1);
 917
 918   // Initialize the index arrays
 919   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 920   lastIndex = new Int_t*[nLayers];
 921   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 922     lastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 923     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 924       lastIndex[iLayer][iModule] =  fLastIndex[iLayer][iModule];
 925       if(iLayer+AliGeomManager::kFirstLayer<=AliGeomManager::kSPD2){
 926         volFit->AddAt(AliGeomManager::LayerToVolUID(iLayer+AliGeomManager::kFirstLayer,iModule),maxntr);
 927         maxntr++;
 928       }
 929     }
 930   }
 931   Int_t found=0;
 932   maxntr=minNtracks+1;
 933   while (maxntr>minNtracks){
 934     maxntr=minNtracks;
 935     for (Int_t iLayer = 0; iLayer <= (AliGeomManager::kSPD2 - AliGeomManager::kFirstLayer); iLayer++) {
 936       for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 937         if(lastIndex[iLayer][iModule]>maxntr){
 938           laymax=iLayer;
 939           modmax=iModule;
 940           maxntr=lastIndex[iLayer][iModule];
 941         }
 942       }
 943     }
 944     if(maxntr>minNtracks){
 945       voluid=AliGeomManager::LayerToVolUID(laymax+AliGeomManager::kFirstLayer,modmax);
 946       sequence->AddAt(voluid,nMod);
 947       lastIndex[laymax][modmax]=0;
 948       nMod++;
 949       volIn->AddAt(voluid,0);
 950     }
 951   }
 952   sequence->Set(nMod);
 953
 954   Int_t ilayer,imod;
 955   for(Int_t iter=0;iter<iterations;iter++){
 956     modAligned=0;
 957     for(Int_t k=0;k<nMod;k++){
 958       TArrayI *volFit3;
 959       voluid=sequence->At(k);
 960       ilayer=AliGeomManager::VolUIDToLayer(voluid,imod);
 961       volIn->AddAt(voluid,0);
 962       found=0;
 963       if(!fitall){
 964         for(Int_t j=0;j<nMod;j++){
 965           if(j==k){
 966             found=1;
 967             continue;
 968           }
 969           else volFit2->AddAt(sequence->At(j),j-found);
 970         }
 971         volFit2->Set(nMod-1);
 972       }
 973       else{
 974         for(Int_t j=0;j<volFit->GetSize();j++){
 975           if(volFit->At(j)!=volIn->At(0))volFit2->AddAt(volFit->At(j),j-found);
 976           else found=1;
 977         }
 978       }
 979
 980       if(volidsSet){
 981         volFit3=IntersectVolArray(volidsSet,volFit2);
 982       }
 983       else volFit3=new TArrayI(*volFit2);
 984
 985
 986       if(AlignVolumesITS(volIn,volFit3,AliGeomManager::kSPD1,AliGeomManager::kSDD1,2))modAligned++;
 987       delete volFit3;
 988       //      if(volidsSet)delete volFit3;
 989     }
 990   }
 991   Int_t noutofsurv=CheckWithSurvey(2.,sequence);
 992   printf("%d modules into the sequence \n %d modules re-aligned \n %d modules moved far away from survey (-> reset) \n",nMod,modAligned,noutofsurv);
 993   delete volFit;
 994   delete volFit2;
 995   delete sequence;
 996
 997   return kTRUE;
 998 }
 999
1000
1001 //__________________________________
1002 Bool_t AliITSRealignTracks::SPDmodulesAlignToSSD(Int_t minNtracks,Int_t iterations){
1003
1004   Int_t volSSD[6]={0,0,0,0,1,1};
1005   TArrayI *volOuter=GetLayersVolUID(volSSD);
1006   TArrayI *voluid=new TArrayI(1);
1007   for(Int_t iter=0;iter<iterations;iter++){
1008     //SPD1
1009     for(Int_t imod=0;imod<AliGeomManager::LayerSize(AliGeomManager::kSPD1);imod++){    if(GetLastIndex(AliGeomManager::kSPD1-AliGeomManager::kFirstLayer,imod)<minNtracks){
1010       printf("Not enough tracks for module: lay %d mod %d \n",1,imod );
1011       continue;
1012     }
1013     voluid->AddAt(AliGeomManager::LayerToVolUID(AliGeomManager::kSPD1,imod),0);
1014     AlignVolumesITS(voluid,volOuter,AliGeomManager::kSSD1,AliGeomManager::kSSD2,2);
1015     }
1016     //SPD2
1017     for(Int_t imod=0;imod<AliGeomManager::LayerSize(AliGeomManager::kSPD2);imod++){
1018       if(GetLastIndex(AliGeomManager::kSPD2-AliGeomManager::kFirstLayer,imod)<minNtracks){
1019         printf("Not enough tracks for module: lay %d mod %d \n",2,imod );
1020         continue;
1021       }
1022       voluid->AddAt(AliGeomManager::LayerToVolUID(AliGeomManager::kSPD2,imod),0);
1023       AlignVolumesITS(voluid,volOuter,AliGeomManager::kSSD1,AliGeomManager::kSSD2,2);
1024     }
1025   }
1026   return kTRUE;
1027 }
1028
1029 //______________________________________________________________________________
1030 Bool_t AliITSRealignTracks::AlignVolumesITS(const TArrayI *volids, const TArrayI *volidsfit,
1031                                      AliGeomManager::ELayerID layerRangeMin,
1032                                      AliGeomManager::ELayerID layerRangeMax,
1033                                      Int_t iterations){
1034
1035   // Align a set of detector volumes.
1036   // Tracks are fitted only within
1037   // the range defined by the user
1038   // (by layerRangeMin and layerRangeMax)
1039   // or within the set of volidsfit
1040   // Repeat the procedure 'iterations' times
1041
1042   Int_t nVolIds = volids->GetSize();
1043   if (nVolIds == 0) {
1044     AliError("Volume IDs array is empty!");
1045     return kFALSE;
1046   }
1047   Bool_t correlated=kFALSE;
1048   Double_t surveycov[21],transl[3],rot[3],survtransl[3],survrot[3];
1049   Double_t frac;
1050
1051   TGeoHMatrix hM;
1052   Double_t phiglob,smearing,rotorig[9],normplanevect[3]={0.,0.,0.},normplanevect2[3]={0.,0.,0.};
1053   Double_t *deltarot;
1054   TMatrixDSym covmatrx(6);
1055   AliAlignObj *modAlign;
1056
1057   // Load only the tracks with at least one
1058   // space point in the set of volume (volids)
1059   BuildIndex();
1060   AliTrackPointArray **points;
1061   Bool_t failed=kFALSE;
1062   Int_t pointsdim=0,skipped=0;
1063   // Start the iterations
1064   while (iterations > 0){
1065     normplanevect2[0]=0.;
1066     normplanevect2[1]=0.;
1067     normplanevect2[2]=0.;
1068     if(fLimitCorr>0.){
1069       ResetCorrModules();
1070       skipped=0;
1071     }
1072     Int_t nArrays = LoadPoints(volids, points,pointsdim);
1073
1074     if (nArrays < fmintracks||nArrays<=0){
1075       failed=kTRUE;
1076       printf("Not enough tracks to try minimization: volUID %d and following in volids \n", volids->At(0));
1077       UnloadPoints(pointsdim, points);
1078       break;
1079     }
1080     frac=1./(Double_t)nArrays;
1081     AliTrackResiduals *minimizer = CreateMinimizer();
1082     minimizer->SetNTracks(nArrays);
1083     minimizer->InitAlignObj();
1084     AliTrackFitter *fitter = CreateFitter();
1085
1086     //Here prepare to set the plane for GetPCArot
1087                                                        //    if(volids->GetSize()==1){//TEMPORARY: to be improved
1088     AliGeomManager::GetOrigRotation(volids->At(0),rotorig);
1089     if((Int_t)AliGeomManager::VolUIDToLayer(volids->At(0))==1){//TEMPORARY: to be improved
1090       normplanevect[0]=-rotorig[1];
1091       normplanevect[1]=-rotorig[4];
1092       normplanevect[2]=0.;
1093     }
1094     else{
1095           normplanevect[0]=rotorig[1];
1096           normplanevect[1]=rotorig[4];
1097           normplanevect[2]=0.;
1098     }
1099
1100     phiglob=TMath::ATan2(normplanevect[1],normplanevect[0]);
1101
1102     modAlign=GetAlignObj(volids->At(0));
1103     modAlign->GetMatrix(hM);
1104     deltarot=hM.GetRotationMatrix();
1105     for(Int_t j=0;j<3;j++){
1106       for(Int_t i=0;i<3;i++){
1107         normplanevect2[j]+=deltarot[j*3+i]*normplanevect[i];
1108       }
1109       // printf("Here the difference: norm1[%d]=%f  norm2[%d]=%f \n",j,normplanevect[j],j,normplanevect2[j]);
1110     }
1111
1112     if(fVarySigmaY){
1113       if(modAlign->GetUniqueID()==0)smearing=fsigmaY;
1114       else{
1115         modAlign->GetCovMatrix(covmatrx);
1116         smearing=5.*5.*(covmatrx(0,0)+covmatrx(1,1)+covmatrx(2,2)+10.*10.*covmatrx(3,3)+10.*10.*covmatrx(4,4)+10.*10.*covmatrx(5,5))/6.;
1117         //This is a sort of average: the trace with the variances of the angles
1118         //weighted with 10 cm divided per 6 and the result multiplied per 25
1119         // (the sqrt would be 5 times a sort of "mean sigma" )
1120         //
1121           }
1122     }
1123     else smearing=fsigmaY;
1124     printf("This is the sigmaY value: %f \n",smearing);
1125     // the plane will be set into the loop on tracks
1126
1127
1128     for (Int_t iArray = 0; iArray < nArrays; iArray++) {
1129       if (!points[iArray]) continue;
1130       points[iArray]->Sort(kTRUE);
1131       fitter->SetTrackPointArray(points[iArray], kFALSE);
1132       // printf("Here normplane vect: %f \n",normplanevect2[1]); //TO BE REPLACED BY      fitter->SetNormPlaneVect(normplanevect2);
1133       if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
1134
1135        if(fLimitCorr>0.){
1136         correlated=kFALSE;
1137         AliTrackPoint p;
1138         Int_t layer,module;
1139         TArrayI *volparray=new TArrayI(points[iArray]->GetNPoints());
1140         for(Int_t point=0;point<points[iArray]->GetNPoints();point++){
1141           points[iArray]->GetPoint(p,point);
1142           volparray->AddAt(p.GetVolumeID(),point);
1143         }
1144         TArrayI *volpArray=ExcludeVolidsFromVolidsArray(volids,volparray);
1145         for(Int_t point=0;point<volpArray->GetSize();point++){
1146           layer=(Int_t)AliGeomManager::VolUIDToLayerSafe(volpArray->At(point),module);
1147           if(fCorrModules[layer-AliGeomManager::kFirstLayer][module]>fLimitCorr){
1148             correlated=kTRUE;
1149             //      printf("volid %d, iarray = %d : skipping %d for Volume: %d \n",volids->At(0),iArray,skipped,volpArray->At(point));
1150             skipped++;
1151             break;
1152           }
1153         }
1154         if(!correlated){
1155           for(Int_t point=0;point<volpArray->GetSize();point++){
1156             layer=(Int_t)AliGeomManager::VolUIDToLayerSafe(volpArray->At(point),module);
1157             //printf("Number of common tracks: %d \n",fCorrModules[layer-AliGeomManager::kFirstLayer][module]);
1158             fCorrModules[layer-AliGeomManager::kFirstLayer][module]+=frac;
1159             delete volparray;
1160             delete volpArray;
1161           }
1162         }
1163         else {
1164           delete volparray;
1165           delete volpArray;
1166           continue;
1167         }
1168        }
1169
1170       AliTrackPointArray *pVolId,*pTrack;
1171       fitter->GetTrackResiduals(pVolId,pTrack);
1172       minimizer->AddTrackPointArrays(pVolId,pTrack);
1173     }
1174
1175     printf("Number of tracks considered: %d \n",nArrays);
1176     frac=(Double_t)skipped/(Double_t)nArrays;
1177     printf("Number of tracks skipped cause of correlation: %d (fraction: %f )\n",skipped,frac);
1178
1179     Int_t ntracks=minimizer->GetNFilledTracks();
1180     frac=(Double_t)ntracks/(Double_t)nArrays;
1181     printf("Number of tracks into the minimizer: %d (fraction: %f )\n",ntracks,frac);
1182     if(ntracks<=fmintracks){
1183       printf("Not enough good tracks found: could not find parameter for volume %d (and following in volids)\n",volids->At(0));
1184       UnloadPoints(pointsdim, points);
1185       failed=kTRUE;
1186       break;
1187     }
1188
1189     failed=(!minimizer->Minimize());
1190
1191     // Update the alignment object(s)
1192     if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
1193       UShort_t volid = (*volids)[iVolId];
1194       if(!failed){
1195         Int_t iModule;
1196         AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
1197         AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
1198
1199         //Check the last minimization is not too large
1200         minimizer->GetAlignObj()->GetPars(transl,rot);
1201         fSurveyObjs[iLayer-AliGeomManager::kFirstLayer][iModule]->GetCovMatrix(surveycov);
1202         fSurveyObjs[iLayer-AliGeomManager::kFirstLayer][iModule]->GetPars(survtransl,survrot);
1203         if(TMath::Sqrt(TMath::Abs(surveycov[0]))*2<TMath::Abs(transl[0]-survtransl[0])||TMath::Sqrt(TMath::Abs(surveycov[2]))*2<TMath::Abs(transl[1]-survtransl[1])||TMath::Sqrt(TMath::Abs(surveycov[5]))*2<TMath::Abs(transl[2]-survtransl[2])||TMath::Sqrt(TMath::Abs(surveycov[9]))*2<TMath::Abs(rot[0]-survrot[0])||TMath::Sqrt(TMath::Abs(surveycov[14]))*2<TMath::Abs(rot[1]-survrot[1])||TMath::Sqrt(TMath::Abs(surveycov[20]))*2<TMath::Abs(rot[2]-survrot[2])){
1204           printf("Results for module %d too large: can't update them \n",volid);
1205           alignObj->SetUniqueID(2);
1206           if(iterations==1){
1207             failed=kTRUE;
1208           }
1209         }
1210         else{
1211           if(fUpdateCov){
1212             *alignObj *= *minimizer->GetAlignObj();
1213             alignObj->SetUniqueID(1);
1214           }
1215           else{
1216             alignObj->GetCovMatrix(covmatrx);
1217             *alignObj *= *minimizer->GetAlignObj();
1218             alignObj->SetCorrMatrix(covmatrx);
1219             alignObj->SetUniqueID(1);
1220           }
1221
1222           /*alignObj->GetPars(transl,rot);
1223
1224           if(TMath::Sqrt(TMath::Abs(surveycov[0]))*20<TMath::Abs(transl[0])||TMath::Sqrt(TMath::Abs(surveycov[2]))*20<TMath::Abs(transl[1])||TMath::Sqrt(TMath::Abs(surveycov[5]))*20<TMath::Abs(transl[2])||TMath::Sqrt(TMath::Abs(surveycov[9]))*20<TMath::Abs(rot[0])||TMath::Sqrt(TMath::Abs(surveycov[14]))*20<TMath::Abs(rot[1])||TMath::Sqrt(TMath::Abs(surveycov[20]))*20<TMath::Abs(rot[2])){
1225           printf("Results for module %d out of Survey: reinitializing it from survey \n",volid);
1226             //    *alignObj = *alignObjSurv;
1227             alignObj->SetPars(0.,0.,0.,0.,0.,0.);
1228             alignObj->SetUniqueID(0);
1229             if(fUpdateCov)alignObj->SetCorrMatrix(surveycov);
1230             if(iterations==1){
1231             failed=kTRUE;
1232             }
1233             }*/
1234         }
1235         if(iterations==1)alignObj->Print("");
1236       }
1237       else {
1238         printf("Minimization failed: cannot update AlignObj for volume: %d \n",volid);
1239       }
1240     }
1241     UnloadPoints(pointsdim,points);
1242     if(failed)break;
1243     minimizer->InitAlignObj();
1244     iterations--;
1245   }
1246
1247   printf("\n \n");
1248
1249   return (!failed);
1250 }
1251
1252
1253
1254 //______________________________________________
1255 Bool_t AliITSRealignTracks::AlignSPDBarrel(Int_t iterations){
1256
1257   Int_t size=0,size2=0;
1258   Int_t layers[6]={1,1,0,0,0,0};
1259   for(Int_t k=1;k<=2;k++){
1260     size+=AliGeomManager::LayerSize(k);
1261   }
1262   for(Int_t k=3;k<=6;k++){
1263     size2+=AliGeomManager::LayerSize(k);
1264     printf("size: %d \n",size2);
1265   }
1266
1267   printf("Aligning SPDBarrel: nmodules: %d \n",size);
1268   printf("Fitting modules: %d \n",size2);
1269
1270   TArrayI *volIDs=GetLayersVolUID(layers);
1271   layers[0]=0;
1272   layers[1]=0;
1273   layers[2]=1;
1274   layers[3]=1;
1275   layers[4]=1;
1276   layers[5]=1;
1277   TArrayI *volIDsFit=GetLayersVolUID(layers);
1278
1279   AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSDD1,AliGeomManager::kTPC1,iterations);
1280
1281   return kTRUE;
1282 }
1283
1284 //______________________
1285 Bool_t AliITSRealignTracks::AlignSPDHalfBarrel(Int_t method,Int_t iterations){
1286
1287   Int_t size=0,size2=0;
1288   Int_t layers[6]={0,0,1,1,1,1};
1289   Int_t sectorsUp[10]={1,1,1,1,1,0,0,0,0,0};
1290   Int_t sectorsDown[10]={0,0,0,0,0,1,1,1,1,1};
1291
1292   TString updownstr;
1293   if(method==0)updownstr="UpNoDown";
1294   else if (method==1)updownstr="DownNoUp";
1295   else if (method==10)updownstr="UpWithDown";
1296   else if (method==11)updownstr="DownWithUp";
1297   else {
1298     AliWarning("Wrong AlignSPDHalfBarrel method selected ");
1299     return kFALSE;
1300   }
1301
1302   for(Int_t i=1;i<=2;i++){
1303     size+=AliGeomManager::LayerSize(i);
1304   }
1305
1306   for(Int_t i=3;i<=6;i++){
1307     size2+=AliGeomManager::LayerSize(i);
1308   }
1309
1310   size=size/2;
1311   if(method==10||method==11)size2+=size;
1312
1313   printf("Aligning  SPDHalfBarrel %s: nmodules: %d \n",updownstr.Data(),size);
1314   printf("Fitting modules: %d \n",size2);
1315   TArrayI *volIDsFit2;
1316   TArrayI *volids = NULL;
1317   TArrayI *volIDsFit=GetLayersVolUID(layers);
1318   if(method==0||method==10)volids=GetSPDSectorsVolids(sectorsUp);
1319   if(method==1||method==11)volids=GetSPDSectorsVolids(sectorsDown);
1320
1321   if(method==10)volIDsFit2=JoinVolArrays(GetSPDSectorsVolids(sectorsDown),volIDsFit);
1322   else if(method==11)volIDsFit2=JoinVolArrays(GetSPDSectorsVolids(sectorsUp),volIDsFit);
1323   else volIDsFit2=volIDsFit;
1324
1325   AlignVolumesITS(volids,volIDsFit2,AliGeomManager::kSPD1,AliGeomManager::kSSD2,iterations);
1326
1327   return kTRUE;
1328 }
1329
1330
1331 //______________________________________________________
1332 Bool_t AliITSRealignTracks::AlignLayer(Int_t layer,Int_t iterations){
1333
1334   Int_t size=0,size2=0;
1335   Int_t layers[6]={0,0,0,0,0,0};
1336   layers[layer-1]=1;
1337   TString layerstr[6]={"SPD1","SPD2","SDD1","SDD2","SSD1","SSD2"};
1338   for(Int_t k=1;k<=6;k++){
1339     if(k!=layer)size2+=AliGeomManager::LayerSize(k);
1340   }
1341   size=AliGeomManager::LayerSize(layer);
1342
1343   printf("Aligning layer %s, nmodules %d ,fitted modules %d \n",layerstr[layer-1].Data(),size,size2);
1344
1345
1346   TArrayI *volIDs=GetLayersVolUID(layers);
1347   layers[0]=1;
1348   layers[1]=1;
1349   layers[2]=1;
1350   layers[3]=1;
1351   layers[4]=1;
1352   layers[5]=1;
1353   layers[layer]=0;
1354   TArrayI *volIDsFit=GetLayersVolUID(layers);
1355
1356   AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSDD1,AliGeomManager::kSSD2,iterations);
1357
1358   return kTRUE;
1359 }
1360
1361 //___________________________________________
1362
1363 Bool_t AliITSRealignTracks::AlignLayersToLayers(Int_t *layer,Int_t iterations){
1364
1365   UShort_t volid;
1366   Int_t size=0,size2=0,j=0,k=0;
1367   Int_t iLayer;
1368   TString layerstr[6]={"SPD1","SPD2","SDD1","SDD2","SSD1","SSD2"};
1369   TString command="",str;
1370   for(Int_t i=1;i<=6;i++){
1371     if(layer[i-1]==1||layer[i-1]==3){
1372       size+=AliGeomManager::LayerSize(i);
1373       command.Append(" ");
1374       command.Append(layerstr[i-1]);
1375     }
1376     if(layer[i-1]==2||layer[i-1]==3){
1377       size2+=AliGeomManager::LayerSize(i);
1378       str.Append(" ");
1379       str.Append(layerstr[i-1]);
1380     }
1381   }
1382
1383   printf("Aligning layers %s To layers %s, nmodules %d ,fitted modules %d \n",command.Data(),str.Data(),size,size2);
1384
1385
1386   TArrayI volIDs(size);
1387   TArrayI volIDsFit(size2);
1388
1389   for (iLayer=(Int_t)AliGeomManager::kSPD1;iLayer<(Int_t)AliGeomManager::kTPC1;iLayer++){
1390     if(layer[iLayer-AliGeomManager::kFirstLayer]==0)continue;
1391     if(layer[iLayer-AliGeomManager::kFirstLayer]==1){
1392       for (Int_t iModule2=0;iModule2<AliGeomManager::LayerSize(iLayer);iModule2++){
1393         volid = AliGeomManager::LayerToVolUID(iLayer,iModule2);
1394         volIDs.AddAt(volid,j);
1395         j++;
1396       }
1397     }
1398     else if(layer[iLayer-AliGeomManager::kFirstLayer]==2){
1399       for (Int_t iModule2=0;iModule2<AliGeomManager::LayerSize(iLayer);iModule2++){
1400         volid = AliGeomManager::LayerToVolUID(iLayer,iModule2);
1401         volIDsFit.AddAt(volid,k);
1402         k++;
1403       }
1404     }
1405     else if(layer[iLayer-AliGeomManager::kFirstLayer]==3){
1406       for (Int_t iModule2=0;iModule2<AliGeomManager::LayerSize(iLayer);iModule2++){
1407         volid = AliGeomManager::LayerToVolUID(iLayer,iModule2);
1408         volIDs.AddAt(volid,j);
1409         j++;
1410         volIDsFit.AddAt(volid,k);
1411         k++;
1412       }
1413     }
1414   }
1415
1416   AlignVolumesITS(&volIDs,&volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,iterations);
1417
1418   return kTRUE;
1419 }
1420
1421 //______________________________________________
1422
1423 Bool_t AliITSRealignTracks::AlignSPDSectorToOuterLayers(Int_t sector,Int_t iterations){
1424
1425
1426   Int_t layers[6]={0,0,1,1,1,1};
1427   Bool_t spd=kFALSE;
1428   Int_t sectorsIN[10]={0,0,0,0,0,0,0,0,0,0};
1429   Int_t sectorsFit[10]={1,1,1,1,1,1,1,1,1,1};
1430
1431   if(sector<0){
1432     sector=-sector;
1433     spd=kTRUE;
1434   }
1435   sectorsIN[sector]=1;
1436   sectorsFit[sector]=0;
1437   TArrayI *volIDs=GetSPDSectorsVolids(sectorsIN);
1438   TArrayI *volIDsFit;
1439   if(spd){
1440     volIDsFit=JoinVolArrays(GetSPDSectorsVolids(sectorsFit),GetLayersVolUID(layers));
1441   }
1442   else volIDsFit=GetLayersVolUID(layers);
1443
1444   printf("Aligning SPD sector %d: nmodules: %d \n",sector,volIDs->GetSize());
1445   printf("Fitting modules: %d \n",volIDsFit->GetSize());
1446
1447   AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSPD2,iterations);
1448
1449   return kTRUE;
1450 }
1451
1452 //______________________________________________
1453 Bool_t AliITSRealignTracks::AlignSPDSectorWithSectors(Int_t sector,Int_t iterations){
1454
1455
1456   Int_t sectorsIN[10]={0,0,0,0,0,0,0,0,0,0};
1457   Int_t sectorsFit[10]={1,1,1,1,1,1,1,1,1,1};
1458
1459   sectorsIN[sector]=1;
1460   sectorsFit[sector]=0;
1461   TArrayI *volIDs=GetSPDSectorsVolids(sectorsIN);
1462   TArrayI *volIDsFit=GetSPDSectorsVolids(sectorsFit);;
1463
1464   printf("Aligning SPD sector %d: nmodules: %d \n",sector,volIDs->GetSize());
1465   printf("Fitting modules: %d \n",volIDsFit->GetSize());
1466
1467
1468   AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSDD1,iterations);
1469
1470   return kTRUE;
1471 }
1472
1473
1474
1475 //___________________________________________________
1476 Bool_t AliITSRealignTracks::AlignSPDSectorsWithSectors(Int_t *sectorsIN,Int_t *sectorsFit,Int_t iterations){
1477
1478   TArrayI *volIDs=GetSPDSectorsVolids(sectorsIN);
1479   TArrayI *volIDsFit=GetSPDSectorsVolids(sectorsFit);;
1480
1481   printf("Aligning SPD sectors: modules: %d \n",volIDs->GetSize());
1482   printf("Fitting modules: %d \n",volIDsFit->GetSize());
1483
1484
1485
1486   return AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSDD1,iterations);;
1487 }
1488
1489 //___________________________________________________
1490 Bool_t AliITSRealignTracks::AlignSPDStaves(Int_t *staves,Int_t *sectorsIN,Int_t *sectorsFit,Int_t iterations){
1491
1492   TArrayI *volIDs=GetSPDStavesVolids(sectorsIN,staves);
1493   TArrayI *volIDsFit=GetSPDSectorsVolids(sectorsFit);
1494
1495   if(volIDs->GetSize()==0){
1496     printf("EMPTY ARRAY !! \n");
1497     return kFALSE;
1498   }
1499   printf("Aligning SPD staves: modules: %d \n",volIDs->GetSize());
1500   printf("Fitting modules: %d \n",volIDsFit->GetSize());
1501
1502   TArrayI *volIDsFit2=ExcludeVolidsFromVolidsArray(volIDs,volIDsFit);
1503   return  AlignVolumesITS(volIDs,volIDsFit2,AliGeomManager::kSPD1,AliGeomManager::kSSD1,iterations);
1504 }
1505
1506
1507 //___________________________________________
1508
1509 Bool_t AliITSRealignTracks::AlignLayerToSPDHalfBarrel(Int_t layer,Int_t updown,Int_t iterations){
1510
1511
1512   Int_t sectorsDown[10]={0,0,0,0,0,1,1,1,1,1};
1513   Int_t sectorsUp[10]={1,1,1,1,1,0,0,0,0,0};
1514   TString layerstr[6]={"SPD1","SPD2","SDD1","SDD2","SSD1","SSD2"};
1515   TArrayI *volIDsFit;
1516   Int_t layers[6]={0,0,0,0,0,0};
1517   layers[layer]=1;
1518   Int_t size=AliGeomManager::LayerSize(layer);
1519   TArrayI *volIDs=GetLayersVolUID(layers);
1520
1521   if(updown==0){
1522     volIDsFit=GetSPDSectorsVolids(sectorsUp);
1523     printf("Aligning layer %s, nmodules %d ,to half barrel Up \n",layerstr[layer-1].Data(),size);
1524   }
1525   else if(updown==1){
1526     volIDsFit=GetSPDSectorsVolids(sectorsDown);
1527     printf("Aligning layer %s, nmodules %d ,to half barrel Down \n",layerstr[layer-1].Data(),size);
1528   }
1529   else {
1530     printf("Wrong Half Barrel selection! \n");
1531     return kFALSE;
1532   }
1533
1534   AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,iterations);
1535
1536   return kTRUE;
1537 }
1538
1539 //___________________________________________
1540
1541 Bool_t AliITSRealignTracks::AlignLayerToSector(Int_t layer,Int_t sector,Int_t iterations){
1542
1543   if(sector>9){
1544     printf("Wrong Sector selection! \n");
1545     return kFALSE;
1546   }
1547   Int_t sectors[10]={0,0,0,0,0,0,0,0,0,0};
1548   sectors[sector]=1;
1549   TString layerstr[6]={"SPD1","SPD2","SDD1","SDD2","SSD1","SSD2"};
1550   TArrayI *volIDsFit;
1551   Int_t layers[6]={0,0,0,0,0,0};
1552   layers[layer]=1;
1553   TArrayI *volIDs=GetLayersVolUID(layers);
1554   Int_t size=AliGeomManager::LayerSize(layer);
1555
1556
1557   volIDsFit=GetSPDSectorsVolids(sectors);
1558   printf("Aligning layer %s, nmodules %d ,to half barrel Up \n",layerstr[layer-1].Data(),size);
1559
1560   AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,iterations);
1561
1562   return kTRUE;
1563 }
1564
1565 //_______________________________________________
1566
1567 Bool_t AliITSRealignTracks::AlignSPDHalfBarrelToHalfBarrel(Int_t updown,Int_t iterations){
1568
1569
1570   Int_t sectorsDown[10]={0,0,0,0,0,1,1,1,1,1};
1571   Int_t sectorsUp[10]={1,1,1,1,1,0,0,0,0,0};
1572
1573   TArrayI *volIDsUp=GetSPDSectorsVolids(sectorsUp);
1574   TArrayI *volIDsDown=GetSPDSectorsVolids(sectorsDown);
1575
1576   if(updown==0){
1577     printf("Aligning SPD HalfBarrel up to half Barrel down : nmodules: %d \n",volIDsUp->GetSize());
1578     printf("Fitting modules: %d \n",volIDsDown->GetSize());
1579     AlignVolumesITS(volIDsUp,volIDsDown,AliGeomManager::kSPD1,AliGeomManager::kSPD2,iterations);
1580   }
1581   else if(updown==1){
1582     printf("Aligning SPD HalfBarrel down to half Barrel Up : nmodules: %d \n",volIDsDown->GetSize());
1583     printf("Fitting modules: %d \n",volIDsUp->GetSize());
1584     AlignVolumesITS(volIDsDown,volIDsUp,AliGeomManager::kSPD1,AliGeomManager::kSPD2,iterations);
1585   }
1586   else {
1587     printf("Wrong Half Barrel selection! \n");
1588     return kFALSE;
1589   }
1590
1591   return kTRUE;
1592 }
1593
1594
1595 //_______________________
1596 Bool_t AliITSRealignTracks::AlignSPDHalfBarrelToSectorRef(Int_t sector,Int_t iterations){
1597
1598   Int_t sectorsIN[10]={0,0,0,0,0,1,1,1,1,1};
1599   Int_t sectorsFit[10]={0,0,0,0,0,0,0,0,0,0};
1600
1601   sectorsFit[sector]=1;
1602
1603   TArrayI *volIDs=GetSPDSectorsVolids(sectorsIN);
1604   TArrayI *volIDsFit=GetSPDSectorsVolids(sectorsFit);
1605
1606   printf("Aligning SPD HalfBarrel to sector 0 %d: nmodules: %d \n",sector,volIDs->GetSize());
1607   printf("Fitting modules: %d \n",volIDsFit->GetSize());
1608
1609
1610   AlignVolumesITS(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSPD2,iterations);
1611
1612
1613   return kTRUE;
1614 }
1615 //_________________________________________
1616 Bool_t AliITSRealignTracks::AlignSPD1SectorRef(Int_t sector,Int_t iterations){
1617
1618   Int_t sectorsIN[10]={0,0,0,0,0,0,0,0,0,0};
1619   Int_t sectorsFit[10]={1,1,1,1,1,1,1,1,1,1};
1620   sectorsIN[sector]=1;
1621   sectorsFit[sector]=0;
1622   TArrayI *volIDs=GetSPDSectorsVolids(sectorsIN);
1623   TArrayI *volIDsFit=GetSPDSectorsVolids(sectorsFit);
1624   Int_t size=volIDs->GetSize();
1625   Int_t size2=volIDsFit->GetSize();
1626   UShort_t volID;
1627   Int_t k=0;
1628
1629   TArrayI *volIDsSPD1=new TArrayI(size-8);
1630   TArrayI *volIDsFit2=new TArrayI(size2+8);
1631
1632   for(Int_t j=0;j<size;j++){
1633     volID=volIDs->At(j);
1634     if(AliGeomManager::VolUIDToLayer(volID)==AliGeomManager::kSPD1){
1635       volIDsSPD1->AddAt(volID,size2+k);
1636       k++;
1637     }
1638     else volIDsFit2->AddAt(volID,j-k);
1639   }
1640
1641
1642   for(Int_t j=0;j<size2;j++){
1643     volID=volIDsFit->At(j);
1644     volIDsFit2->AddAt(volID,size-k+j);
1645   }
1646
1647   printf("Aligning SPD Sector %d: nmodules: %d \n",sector,volIDsSPD1->GetSize());
1648   printf("Fitting modules: %d \n",volIDsFit2->GetSize());
1649
1650   AlignVolumesITS(volIDsSPD1,volIDsFit2,AliGeomManager::kSPD1,AliGeomManager::kSPD2,iterations);
1651
1652   return kTRUE;
1653 }
1654
1655 //_____________________________________________
1656
1657 AliAlignObjParams* AliITSRealignTracks::MediateAlignObj(TArrayI *volIDs,Int_t lastVolid){
1658   UShort_t volid;
1659
1660   TGeoHMatrix hm;
1661   Double_t *rot,*transl;
1662   Double_t rotSum[9],translSum[3]={0.,0.,0.};
1663   for(Int_t k=0;k<8;k++)rotSum[k]=0.;
1664
1665
1666   for(Int_t ivol=0;ivol<lastVolid;ivol++){
1667     volid=volIDs->At(ivol);
1668
1669     GetAlignObj(volIDs->At(ivol))->GetMatrix(hm);
1670
1671     rot=hm.GetRotationMatrix();
1672     transl=hm.GetTranslation();
1673
1674     for(Int_t j=0;j<9;j++)rotSum[j]+=rot[j];
1675     for(Int_t jt=0;jt<3;jt++)translSum[jt]+=transl[jt];
1676   }
1677   if(lastVolid!=0){
1678     for(Int_t j=0;j<9;j++)rotSum[j]=rotSum[j]/lastVolid;
1679     for(Int_t jt=0;jt<3;jt++)translSum[jt]=translSum[jt]/lastVolid;
1680   }
1681   else printf("Try to mediate results for zero modules \n");
1682
1683   hm.SetRotation(rotSum);
1684   hm.SetTranslation(translSum);
1685
1686
1687
1688   AliAlignObjParams *alignObj=new AliAlignObjParams("average", 0,hm, kTRUE);
1689   return alignObj;
1690
1691 }
1692
1693
1694 //________________________________________________
1695 TArrayI* AliITSRealignTracks::GetSPDStavesVolids(Int_t *sectors,Int_t* staves){
1696
1697
1698   // This method gets the volID Array for the chosen staves into the
1699   // chosen sectors. You have to pass an array (10 dim) with a 1 for each
1700   // selected sector and an array (6 dim) with a 1 for each chosen stave.
1701   // The staves are numbered in this way: 0,1 for SPD1 and 2,3,4,5 for SPD2
1702   // i.e. sectors[10] = {1,1,0,0,0,0,0,0,1,0} -> Sector 0, 1, 9 selected.
1703   // staves[6]={0,1,1,0,0,1} -> Staves 1 on SPD1 and 0 and 3 on SPD2 selected
1704
1705   Int_t nSect=0,nStaves=0;
1706   Int_t last=0;
1707
1708
1709   for(Int_t co=0;co<10;co++){ //counts the number of sectors chosen
1710     if(sectors[co]==1) nSect++;
1711   }
1712
1713   for(Int_t co=0;co<6;co++){ //counts the number of sectors chosen
1714     if(staves[co]==1) nStaves++;
1715   }
1716
1717   if(nSect<1||nStaves<1){ //if no sector chosen -> exit
1718     Printf("Error! No Sector/s or staves Selected!");
1719     return 0x0;
1720   }
1721
1722   TArrayI *volIDs = new TArrayI(nSect*nStaves*4);
1723   TString stave="/Stave",str,symn,laystr;
1724
1725   TArrayI *sectvol=GetSPDSectorsVolids(sectors);
1726   //SPD1
1727   laystr="SPD0";
1728   for(Int_t k=0;k<2;k++){
1729     if(staves[k]==1){
1730       str=stave;
1731       str+=k;
1732       for(Int_t i=0;i<sectvol->GetSize();i++){
1733         symn=AliGeomManager::SymName(sectvol->At(i));
1734         if(symn.Contains(str)&&symn.Contains(laystr)){
1735           //      printf("Adding: %s \n",symn.Data());
1736           volIDs->AddAt(sectvol->At(i),last);
1737           last++;
1738         }
1739       }
1740     }
1741   }
1742   //SPD1
1743   laystr="SPD1";
1744   for(Int_t k=2;k<6;k++){
1745     if(staves[k]==1){
1746       str=stave;
1747       str+=k-2;
1748       for(Int_t i=0;i<sectvol->GetSize();i++){
1749         symn=AliGeomManager::SymName(sectvol->At(i));
1750         if(symn.Contains(str)&&symn.Contains(laystr)){
1751           volIDs->AddAt(sectvol->At(i),last);
1752           printf("Adding: %s \n",symn.Data());
1753           last++;
1754         }
1755       }
1756     }
1757   }
1758
1759   volIDs->Set(last);
1760   return volIDs;
1761 }
1762
1763 //________________________________________________
1764 TArrayI* AliITSRealignTracks::GetSPDSectorsVolids(Int_t *sectors)
1765 {
1766   //
1767   // This method gets the volID Array for the chosen sectors.
1768   // You have to pass an array with a 1 for each selected sector.
1769   // i.e. sectors[10] = {1,1,0,0,0,0,0,0,1,0} -> Sector 0, 1, 9 selected.
1770   //
1771
1772   Int_t nSect=0;
1773   Int_t iModule=0;
1774
1775
1776   for(Int_t co=0;co<10;co++){ //counts the number of sectors chosen
1777     if(sectors[co]==1) nSect++;
1778   }
1779
1780   if(nSect<1){ //if no sector chosen -> exit
1781     Printf("Error! No Sector/s Selected!");
1782     return 0x0;
1783   }
1784
1785   TArrayI *volIDs = new TArrayI(nSect*24);
1786
1787     if(sectors[0]==1){ //--->cSect = 0 <---
1788       volIDs->AddAt(2048,iModule); iModule++;
1789       volIDs->AddAt(2049,iModule); iModule++;
1790       volIDs->AddAt(2050,iModule); iModule++;
1791       volIDs->AddAt(2051,iModule); iModule++;
1792       volIDs->AddAt(2052,iModule); iModule++;
1793       volIDs->AddAt(2053,iModule); iModule++;
1794       volIDs->AddAt(2054,iModule); iModule++;
1795       volIDs->AddAt(2055,iModule); iModule++;
1796       volIDs->AddAt(4096,iModule); iModule++;
1797       volIDs->AddAt(4097,iModule); iModule++;
1798       volIDs->AddAt(4098,iModule); iModule++;
1799       volIDs->AddAt(4099,iModule); iModule++;
1800       volIDs->AddAt(4100,iModule); iModule++;
1801       volIDs->AddAt(4101,iModule); iModule++;
1802       volIDs->AddAt(4102,iModule); iModule++;
1803       volIDs->AddAt(4103,iModule); iModule++;
1804       volIDs->AddAt(4104,iModule); iModule++;
1805       volIDs->AddAt(4105,iModule); iModule++;
1806       volIDs->AddAt(4106,iModule); iModule++;
1807       volIDs->AddAt(4107,iModule); iModule++;
1808       volIDs->AddAt(4108,iModule); iModule++;
1809       volIDs->AddAt(4109,iModule); iModule++;
1810       volIDs->AddAt(4110,iModule); iModule++;
1811       volIDs->AddAt(4111,iModule); iModule++;
1812     }
1813     if(sectors[1]==1){ //--->cSect = 1 <//---
1814       volIDs->AddAt(2056,iModule); iModule++;
1815       volIDs->AddAt(2057,iModule); iModule++;
1816       volIDs->AddAt(2058,iModule); iModule++;
1817       volIDs->AddAt(2059,iModule); iModule++;
1818       volIDs->AddAt(2060,iModule); iModule++;
1819       volIDs->AddAt(2061,iModule); iModule++;
1820       volIDs->AddAt(2062,iModule); iModule++;
1821       volIDs->AddAt(2063,iModule); iModule++;
1822       volIDs->AddAt(4112,iModule); iModule++;
1823       volIDs->AddAt(4113,iModule); iModule++;
1824       volIDs->AddAt(4114,iModule); iModule++;
1825       volIDs->AddAt(4115,iModule); iModule++;
1826       volIDs->AddAt(4116,iModule); iModule++;
1827       volIDs->AddAt(4117,iModule); iModule++;
1828       volIDs->AddAt(4118,iModule); iModule++;
1829       volIDs->AddAt(4119,iModule); iModule++;
1830       volIDs->AddAt(4120,iModule); iModule++;
1831       volIDs->AddAt(4121,iModule); iModule++;
1832       volIDs->AddAt(4122,iModule); iModule++;
1833       volIDs->AddAt(4123,iModule); iModule++;
1834       volIDs->AddAt(4124,iModule); iModule++;
1835       volIDs->AddAt(4125,iModule); iModule++;
1836       volIDs->AddAt(4126,iModule); iModule++;
1837       volIDs->AddAt(4127,iModule); iModule++;
1838     }
1839     if(sectors[2]==1){//--->cSect = 2 <//---
1840       volIDs->AddAt(2064,iModule); iModule++;
1841       volIDs->AddAt(2065,iModule); iModule++;
1842       volIDs->AddAt(2066,iModule); iModule++;
1843       volIDs->AddAt(2067,iModule); iModule++;
1844       volIDs->AddAt(2068,iModule); iModule++;
1845       volIDs->AddAt(2069,iModule); iModule++;
1846       volIDs->AddAt(2070,iModule); iModule++;
1847       volIDs->AddAt(2071,iModule); iModule++;
1848       volIDs->AddAt(4128,iModule); iModule++;
1849       volIDs->AddAt(4129,iModule); iModule++;
1850       volIDs->AddAt(4130,iModule); iModule++;
1851       volIDs->AddAt(4131,iModule); iModule++;
1852       volIDs->AddAt(4132,iModule); iModule++;
1853       volIDs->AddAt(4133,iModule); iModule++;
1854       volIDs->AddAt(4134,iModule); iModule++;
1855       volIDs->AddAt(4135,iModule); iModule++;
1856       volIDs->AddAt(4136,iModule); iModule++;
1857       volIDs->AddAt(4137,iModule); iModule++;
1858       volIDs->AddAt(4138,iModule); iModule++;
1859       volIDs->AddAt(4139,iModule); iModule++;
1860       volIDs->AddAt(4140,iModule); iModule++;
1861       volIDs->AddAt(4141,iModule); iModule++;
1862       volIDs->AddAt(4142,iModule); iModule++;
1863       volIDs->AddAt(4143,iModule); iModule++;
1864     }
1865     if(sectors[3]==1){//--->cSect = 3 <//---
1866       volIDs->AddAt(2072,iModule); iModule++;
1867       volIDs->AddAt(2073,iModule); iModule++;
1868       volIDs->AddAt(2074,iModule); iModule++;
1869       volIDs->AddAt(2075,iModule); iModule++;
1870       volIDs->AddAt(2076,iModule); iModule++;
1871       volIDs->AddAt(2077,iModule); iModule++;
1872       volIDs->AddAt(2078,iModule); iModule++;
1873       volIDs->AddAt(2079,iModule); iModule++;
1874       volIDs->AddAt(4144,iModule); iModule++;
1875       volIDs->AddAt(4145,iModule); iModule++;
1876       volIDs->AddAt(4146,iModule); iModule++;
1877       volIDs->AddAt(4147,iModule); iModule++;
1878       volIDs->AddAt(4148,iModule); iModule++;
1879       volIDs->AddAt(4149,iModule); iModule++;
1880       volIDs->AddAt(4150,iModule); iModule++;
1881       volIDs->AddAt(4151,iModule); iModule++;
1882       volIDs->AddAt(4152,iModule); iModule++;
1883       volIDs->AddAt(4153,iModule); iModule++;
1884       volIDs->AddAt(4154,iModule); iModule++;
1885       volIDs->AddAt(4155,iModule); iModule++;
1886       volIDs->AddAt(4156,iModule); iModule++;
1887       volIDs->AddAt(4157,iModule); iModule++;
1888       volIDs->AddAt(4158,iModule); iModule++;
1889       volIDs->AddAt(4159,iModule); iModule++;
1890     }
1891     if(sectors[4]==1){//--->cSect = 4 <//---
1892       volIDs->AddAt(2080,iModule); iModule++;
1893       volIDs->AddAt(2081,iModule); iModule++;
1894       volIDs->AddAt(2082,iModule); iModule++;
1895       volIDs->AddAt(2083,iModule); iModule++;
1896       volIDs->AddAt(2084,iModule); iModule++;
1897       volIDs->AddAt(2085,iModule); iModule++;
1898       volIDs->AddAt(2086,iModule); iModule++;
1899       volIDs->AddAt(2087,iModule); iModule++;
1900       volIDs->AddAt(4160,iModule); iModule++;
1901       volIDs->AddAt(4161,iModule); iModule++;
1902       volIDs->AddAt(4162,iModule); iModule++;
1903       volIDs->AddAt(4163,iModule); iModule++;
1904       volIDs->AddAt(4164,iModule); iModule++;
1905       volIDs->AddAt(4165,iModule); iModule++;
1906       volIDs->AddAt(4166,iModule); iModule++;
1907       volIDs->AddAt(4167,iModule); iModule++;
1908       volIDs->AddAt(4168,iModule); iModule++;
1909       volIDs->AddAt(4169,iModule); iModule++;
1910       volIDs->AddAt(4170,iModule); iModule++;
1911       volIDs->AddAt(4171,iModule); iModule++;
1912       volIDs->AddAt(4172,iModule); iModule++;
1913       volIDs->AddAt(4173,iModule); iModule++;
1914       volIDs->AddAt(4174,iModule); iModule++;
1915       volIDs->AddAt(4175,iModule); iModule++;
1916     }
1917     if(sectors[5]==1){//--->cSect = 5 <//---
1918       volIDs->AddAt(2088,iModule); iModule++;
1919       volIDs->AddAt(2089,iModule); iModule++;
1920       volIDs->AddAt(2090,iModule); iModule++;
1921       volIDs->AddAt(2091,iModule); iModule++;
1922       volIDs->AddAt(2092,iModule); iModule++;
1923       volIDs->AddAt(2093,iModule); iModule++;
1924       volIDs->AddAt(2094,iModule); iModule++;
1925       volIDs->AddAt(2095,iModule); iModule++;
1926       volIDs->AddAt(4176,iModule); iModule++;
1927       volIDs->AddAt(4177,iModule); iModule++;
1928       volIDs->AddAt(4178,iModule); iModule++;
1929       volIDs->AddAt(4179,iModule); iModule++;
1930       volIDs->AddAt(4180,iModule); iModule++;
1931       volIDs->AddAt(4181,iModule); iModule++;
1932       volIDs->AddAt(4182,iModule); iModule++;
1933       volIDs->AddAt(4183,iModule); iModule++;
1934       volIDs->AddAt(4184,iModule); iModule++;
1935       volIDs->AddAt(4185,iModule); iModule++;
1936       volIDs->AddAt(4186,iModule); iModule++;
1937       volIDs->AddAt(4187,iModule); iModule++;
1938       volIDs->AddAt(4188,iModule); iModule++;
1939       volIDs->AddAt(4189,iModule); iModule++;
1940       volIDs->AddAt(4190,iModule); iModule++;
1941       volIDs->AddAt(4191,iModule); iModule++;
1942     }
1943     if(sectors[6]==1){//--->cSect = 6 <//---
1944       volIDs->AddAt(2096,iModule); iModule++;
1945       volIDs->AddAt(2097,iModule); iModule++;
1946       volIDs->AddAt(2098,iModule); iModule++;
1947       volIDs->AddAt(2099,iModule); iModule++;
1948       volIDs->AddAt(2100,iModule); iModule++;
1949       volIDs->AddAt(2101,iModule); iModule++;
1950       volIDs->AddAt(2102,iModule); iModule++;
1951       volIDs->AddAt(2103,iModule); iModule++;
1952       volIDs->AddAt(4192,iModule); iModule++;
1953       volIDs->AddAt(4193,iModule); iModule++;
1954       volIDs->AddAt(4194,iModule); iModule++;
1955       volIDs->AddAt(4195,iModule); iModule++;
1956       volIDs->AddAt(4196,iModule); iModule++;
1957       volIDs->AddAt(4197,iModule); iModule++;
1958       volIDs->AddAt(4198,iModule); iModule++;
1959       volIDs->AddAt(4199,iModule); iModule++;
1960       volIDs->AddAt(4200,iModule); iModule++;
1961       volIDs->AddAt(4201,iModule); iModule++;
1962       volIDs->AddAt(4202,iModule); iModule++;
1963       volIDs->AddAt(4203,iModule); iModule++;
1964       volIDs->AddAt(4204,iModule); iModule++;
1965       volIDs->AddAt(4205,iModule); iModule++;
1966       volIDs->AddAt(4206,iModule); iModule++;
1967       volIDs->AddAt(4207,iModule); iModule++;
1968     }
1969      if(sectors[7]==1){ //--->cSect = 7 <//---
1970        volIDs->AddAt(2104,iModule); iModule++;
1971        volIDs->AddAt(2105,iModule); iModule++;
1972        volIDs->AddAt(2106,iModule); iModule++;
1973        volIDs->AddAt(2107,iModule); iModule++;
1974        volIDs->AddAt(2108,iModule); iModule++;
1975        volIDs->AddAt(2109,iModule); iModule++;
1976        volIDs->AddAt(2110,iModule); iModule++;
1977        volIDs->AddAt(2111,iModule); iModule++;
1978        volIDs->AddAt(4208,iModule); iModule++;
1979        volIDs->AddAt(4209,iModule); iModule++;
1980        volIDs->AddAt(4210,iModule); iModule++;
1981        volIDs->AddAt(4211,iModule); iModule++;
1982        volIDs->AddAt(4212,iModule); iModule++;
1983        volIDs->AddAt(4213,iModule); iModule++;
1984        volIDs->AddAt(4214,iModule); iModule++;
1985        volIDs->AddAt(4215,iModule); iModule++;
1986        volIDs->AddAt(4216,iModule); iModule++;
1987        volIDs->AddAt(4217,iModule); iModule++;
1988        volIDs->AddAt(4218,iModule); iModule++;
1989        volIDs->AddAt(4219,iModule); iModule++;
1990        volIDs->AddAt(4220,iModule); iModule++;
1991        volIDs->AddAt(4221,iModule); iModule++;
1992        volIDs->AddAt(4222,iModule); iModule++;
1993        volIDs->AddAt(4223,iModule); iModule++;
1994      }
1995      if(sectors[8]==1){//--->cSect = 8 <//---
1996        volIDs->AddAt(2112,iModule); iModule++;
1997        volIDs->AddAt(2113,iModule); iModule++;
1998        volIDs->AddAt(2114,iModule); iModule++;
1999        volIDs->AddAt(2115,iModule); iModule++;
2000        volIDs->AddAt(2116,iModule); iModule++;
2001        volIDs->AddAt(2117,iModule); iModule++;
2002        volIDs->AddAt(2118,iModule); iModule++;
2003        volIDs->AddAt(2119,iModule); iModule++;
2004        volIDs->AddAt(4224,iModule); iModule++;
2005        volIDs->AddAt(4225,iModule); iModule++;
2006        volIDs->AddAt(4226,iModule); iModule++;
2007        volIDs->AddAt(4227,iModule); iModule++;
2008        volIDs->AddAt(4228,iModule); iModule++;
2009        volIDs->AddAt(4229,iModule); iModule++;
2010        volIDs->AddAt(4230,iModule); iModule++;
2011        volIDs->AddAt(4231,iModule); iModule++;
2012        volIDs->AddAt(4232,iModule); iModule++;
2013        volIDs->AddAt(4233,iModule); iModule++;
2014        volIDs->AddAt(4234,iModule); iModule++;
2015        volIDs->AddAt(4235,iModule); iModule++;
2016        volIDs->AddAt(4236,iModule); iModule++;
2017        volIDs->AddAt(4237,iModule); iModule++;
2018        volIDs->AddAt(4238,iModule); iModule++;
2019        volIDs->AddAt(4239,iModule); iModule++;
2020      }
2021      if(sectors[9]==1){//--->cSect = 9 <//---
2022        volIDs->AddAt(2120,iModule); iModule++;
2023        volIDs->AddAt(2121,iModule); iModule++;
2024        volIDs->AddAt(2122,iModule); iModule++;
2025        volIDs->AddAt(2123,iModule); iModule++;
2026        volIDs->AddAt(2124,iModule); iModule++;
2027        volIDs->AddAt(2125,iModule); iModule++;
2028        volIDs->AddAt(2126,iModule); iModule++;
2029        volIDs->AddAt(2127,iModule); iModule++;
2030        volIDs->AddAt(4240,iModule); iModule++;
2031        volIDs->AddAt(4241,iModule); iModule++;
2032        volIDs->AddAt(4242,iModule); iModule++;
2033        volIDs->AddAt(4243,iModule); iModule++;
2034        volIDs->AddAt(4244,iModule); iModule++;
2035        volIDs->AddAt(4245,iModule); iModule++;
2036        volIDs->AddAt(4246,iModule); iModule++;
2037        volIDs->AddAt(4247,iModule); iModule++;
2038        volIDs->AddAt(4248,iModule); iModule++;
2039        volIDs->AddAt(4249,iModule); iModule++;
2040        volIDs->AddAt(4250,iModule); iModule++;
2041        volIDs->AddAt(4251,iModule); iModule++;
2042        volIDs->AddAt(4252,iModule); iModule++;
2043        volIDs->AddAt(4253,iModule); iModule++;
2044        volIDs->AddAt(4254,iModule); iModule++;
2045        volIDs->AddAt(4255,iModule); iModule++;
2046      }
2047
2048   return volIDs;
2049 }
2050
2051 //___________________________________
2052 TArrayI* AliITSRealignTracks::GetLayersVolUID(Int_t *layer){
2053
2054   TArrayI *out=new TArrayI(2198);
2055   Int_t last=0;
2056   UShort_t voluid;
2057   for(Int_t i=0;i<6;i++){
2058     if(layer[i]==1){
2059       for(Int_t mod=0;mod<AliGeomManager::LayerSize(i+AliGeomManager::kFirstLayer);mod++){
2060         voluid=AliGeomManager::LayerToVolUID(i+AliGeomManager::kFirstLayer,mod);
2061         out->AddAt(voluid,last);
2062         //      printf("voluid %d at position %d \n",out->At(last),last);
2063         last++;
2064       }
2065     }
2066   }
2067   out->Set(last);
2068   return out;
2069 }
2070
2071 //_________________
2072 TArrayI* AliITSRealignTracks::SelectLayerInVolids(const TArrayI *volidsIN,AliGeomManager::ELayerID layer){
2073
2074   Int_t size=volidsIN->GetSize();
2075   Int_t count=0;
2076   for(Int_t j=0;j<size;j++){
2077     if(AliGeomManager::VolUIDToLayer(volidsIN->At(j))==layer)count++;
2078   }
2079   TArrayI *volidsOUT=new TArrayI(count);
2080   count=0;
2081   for(Int_t j=0;j<size;j++){
2082     if(AliGeomManager::VolUIDToLayer(volidsIN->At(j))==layer){
2083       volidsOUT->AddAt(volidsIN->At(j),count);
2084       count++;
2085     }
2086   }
2087   return volidsOUT;
2088 }
2089
2090 //______________________________________________
2091
2092 TArrayI* AliITSRealignTracks::IntersectVolArray(const TArrayI *vol1,const TArrayI *vol2){
2093
2094   Int_t size1=vol1->GetSize();
2095   Int_t size2=vol2->GetSize();
2096   Int_t last=0,volid;
2097   Bool_t found;
2098   TArrayI *volidOut=new TArrayI(size1+size2);
2099
2100   for(Int_t k=0;k<size1;k++){
2101     found=kFALSE;
2102     volid=vol1->At(k);
2103     for(Int_t j=0;j<size2;j++){
2104       if(vol2->At(j)==volid)found=kTRUE;
2105     }
2106     if(found){
2107       volidOut->AddAt(volid,last);
2108       last++;
2109     }
2110   }
2111   volidOut->Set(last);
2112   return volidOut;
2113 }
2114 //_________________________________________
2115
2116 TArrayI* AliITSRealignTracks::JoinVolArrays(const TArrayI *vol1,const TArrayI *vol2){
2117   //!BE CAREFUL: If an index is repeated in vol1 or vol2 will be repeated also in the final array
2118
2119   Int_t size1=vol1->GetSize();
2120   Int_t size2=vol2->GetSize();
2121   Int_t count=0;
2122   UShort_t volid;
2123   Bool_t found;
2124   TArrayI *volidOut=new TArrayI(size1+size2);
2125
2126   for(Int_t k=0;k<size1;k++){
2127     volid=vol1->At(k);
2128     volidOut->AddAt(volid,k);
2129   }
2130
2131   for(Int_t k=0;k<size1;k++){
2132     found=kFALSE;
2133     volid=vol1->At(k);
2134     for(Int_t j=0;j<size2;j++){
2135       if(vol2->At(j)==volid)found=kTRUE;
2136     }
2137     if(!found){
2138       volidOut->AddAt(volid,size1+count);
2139       count++;
2140     }
2141   }
2142   volidOut->Set(size1+count);
2143   return volidOut;
2144 }
2145
2146 //______________________________________
2147
2148 TArrayI* AliITSRealignTracks::ExcludeVolidsFromVolidsArray(const TArrayI *volidsToExclude,const TArrayI *volStart){
2149
2150   Int_t size1=volidsToExclude->GetSize();
2151   Int_t size2=volStart->GetSize();
2152   Int_t last=0;
2153   UShort_t volid;
2154   Bool_t found;
2155   TArrayI *volidOut=new TArrayI(size2);
2156
2157   for(Int_t k=0;k<size2;k++){
2158     found=kFALSE;
2159     volid=volStart->At(k);
2160     for(Int_t j=0;j<size1;j++){
2161       if(volidsToExclude->At(j)==volid){
2162         found=kTRUE;
2163         break;
2164       }
2165     }
2166     if(!found){
2167       volidOut->AddAt(volid,last);
2168       last++;
2169     }
2170   }
2171   volidOut->Set(last);
2172   return volidOut;
2173 }
2174
2175
2176 //________________________________________
2177
2178 TArrayI* AliITSRealignTracks::GetLayerVolumes(Int_t *layer){
2179
2180   TArrayI *out=new TArrayI(2198);
2181   Int_t last=0;
2182   UShort_t voluid;
2183   for(Int_t i=0;i<6;i++){
2184     if(layer[i]==1){
2185       for(Int_t mod=0;mod<AliGeomManager::LayerSize(i+AliGeomManager::kFirstLayer);mod++){
2186         voluid=AliGeomManager::LayerToVolUID(i+AliGeomManager::kFirstLayer,mod);
2187         out->AddAt(voluid,last);
2188         //      printf("voluid %d at position %d \n",out->At(last),last);
2189         last++;
2190       }
2191     }
2192   }
2193   out->Set(last);
2194   return out;
2195 }
2196
2197
2198
2199 //______________________________
2200 TArrayI* AliITSRealignTracks::GetAlignedVolumes(char *filename){
2201   if(gSystem->AccessPathName(filename)){
2202     printf("Wrong Realignment file name \n");
2203     return 0x0;
2204   }
2205   TFile *f=TFile::Open(filename,"READ");
2206   TClonesArray *array=(TClonesArray*)f->Get("ITSAlignObjs");
2207   AliAlignObjParams *a;
2208   Int_t last=0;
2209   TArrayI *volidOut=new TArrayI(2200);
2210   for(Int_t j=0;j<array->GetSize();j++){
2211     a=(AliAlignObjParams*)array->At(j);
2212     if(a->GetUniqueID()==0)continue;
2213
2214     else {
2215       volidOut->AddAt(a->GetVolUID(),last);
2216       last++;
2217     }
2218   }
2219   volidOut->Set(last);
2220   f->Close();
2221   return volidOut;
2222 }
2223
2224
2225 //________________________________________
2226 void AliITSRealignTracks::SetDraw(Bool_t draw,Bool_t refresh){
2227
2228   if(refresh){
2229     // WriteHists();
2230     if(fAlignDrawObjs)DeleteDrawHists();
2231     InitDrawHists();
2232   }
2233   fDraw=draw;
2234   return;
2235 }
2236
2237 void AliITSRealignTracks::DeleteDrawHists(){
2238   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
2239     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
2240       delete fAlignDrawObjs[iLayer][iModule];
2241     }
2242     if(fAlignDrawObjs[iLayer])delete [] fAlignDrawObjs[iLayer];
2243   }
2244
2245   delete [] fAlignDrawObjs;
2246   fAlignDrawObjs = 0;
2247
2248
2249   delete fCanvPar;
2250   delete fCanvGr;
2251   delete fgrIterMeanX;
2252   delete fgrIterRMSX;
2253   delete fgrIterMeanY;
2254   delete fgrIterRMSY;
2255   delete fgrIterMeanZ;
2256   delete fgrIterRMSZ;
2257   delete fgrIterMeanPsi;
2258   delete fgrIterRMSPsi;
2259   delete fgrIterMeanTheta;
2260   delete fgrIterRMSTheta;
2261   delete fgrIterMeanPhi;
2262   delete fgrIterRMSPhi;
2263
2264 }
2265
2266 void AliITSRealignTracks::InitDrawHists(){
2267   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
2268   fAlignDrawObjs = new AliAlignObj**[nLayers];
2269   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
2270     fAlignDrawObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
2271     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
2272       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
2273       fAlignDrawObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
2274       fAlignDrawObjs[iLayer][iModule]->SetUniqueID(1);
2275     }
2276   }
2277
2278
2279   TH1F *hX=new TH1F("hX","hX",1000,-10000.,10000.);
2280   TH1F *hY=new TH1F("hY","hY",1000,-10000.,10000.);
2281   TH1F *hZ=new TH1F("hZ","hZ",1000,-10000.,10000.);
2282   TH1F *hPsi=new TH1F("hPsi","hPsi",1000,-5000.,5000.);
2283   TH1F *hTheta=new TH1F("hTheta","hTheta",1000,-5000.,5000.);
2284   TH1F *hPhi=new TH1F("hPhi","hPhi",1000,-5000.,5000.);
2285
2286   fCanvPar=new TCanvas("fCanvPar","Parameters trend during iterations: Convergence \n");
2287   fCanvPar->Divide(3,2);
2288   fCanvPar->cd(1);
2289   hX->Draw();
2290   hX->SetXTitle("#mum");
2291   fCanvPar->cd(2);
2292   hY->Draw();
2293   hY->SetXTitle("#mum");
2294   fCanvPar->cd(3);
2295   hZ->SetXTitle("#mum");
2296   hZ->Draw();
2297   fCanvPar->cd(4);
2298   hPsi->SetXTitle("mdeg");
2299   hPsi->Draw();
2300   fCanvPar->cd(5);
2301   hTheta->SetXTitle("mdeg");
2302   hTheta->Draw();
2303   fCanvPar->cd(6);
2304   hPhi->SetXTitle("mdeg");
2305   hPhi->Draw();
2306   fCanvPar->Update();
2307
2308
2309   fCanvGr=new TCanvas("fCanvGr","Parameters trend during iterations: Convergence \n");
2310   fCanvGr->Divide(3,2);
2311
2312   fCanvGr->cd(1);
2313   fgrIterMeanX=new TGraph(1);
2314   fgrIterRMSX=new TGraph(1);
2315   fgrIterRMSX->GetYaxis()->SetRangeUser(-1000.,1000.);
2316   fgrIterRMSX->SetName("fgrIterRMSX");
2317   fgrIterRMSX->SetLineColor(2);
2318   fgrIterMeanX->SetName("fgrIterMeanX");
2319   fgrIterMeanX->SetTitle("Convergence of #deltaX \n");
2320   fgrIterMeanX->GetXaxis()->SetTitle("#mum");
2321   fgrIterRMSX->Draw("acp");
2322   fgrIterMeanX->Draw("cp");
2323
2324   fCanvGr->cd(2);
2325   fgrIterMeanY=new TGraph(1);
2326   fgrIterRMSY=new TGraph(1);
2327   fgrIterRMSY->GetYaxis()->SetRangeUser(-1000.,1000.);
2328   fgrIterRMSY->SetName("fgrIterRMSY");
2329   fgrIterRMSY->SetLineColor(2);
2330   fgrIterMeanY->SetName("fgrIterMeanY");
2331   fgrIterMeanY->SetTitle("Convergence of #deltaY \n");
2332   fgrIterMeanY->GetXaxis()->SetTitle("#mum");
2333   fgrIterRMSY->Draw("acp");
2334   fgrIterMeanY->Draw("cp");
2335
2336   fCanvGr->cd(3);
2337   fgrIterMeanZ=new TGraph(1);
2338   fgrIterRMSZ=new TGraph(1);
2339   fgrIterRMSZ->GetYaxis()->SetRangeUser(-1000.,1000.);
2340   fgrIterRMSZ->SetName("fgrIterRMSZ");
2341   fgrIterRMSZ->SetLineColor(2);
2342   fgrIterMeanZ->SetName("fgrIterMeanZ");
2343   fgrIterMeanZ->SetTitle("Convergence of #deltaZ \n");
2344   fgrIterMeanZ->GetXaxis()->SetTitle("#mum");
2345   fgrIterRMSZ->Draw("acp");
2346   fgrIterMeanZ->Draw("cp");
2347
2348   fCanvGr->cd(4);
2349   fgrIterMeanPsi=new TGraph(1);
2350   fgrIterRMSPsi=new TGraph(1);
2351   fgrIterRMSPsi->GetYaxis()->SetRangeUser(-1000.,1000.);
2352   fgrIterRMSPsi->SetName("fgrIterRMSPsi");
2353   fgrIterRMSPsi->SetLineColor(2);
2354   fgrIterMeanPsi->SetName("fgrIterMeanPsi");
2355   fgrIterMeanPsi->SetTitle("Convergence of #deltaPsi \n");
2356   fgrIterMeanPsi->GetXaxis()->SetTitle("mdeg");
2357   fgrIterRMSPsi->Draw("acp");
2358   fgrIterMeanPsi->Draw("cp");
2359
2360   fCanvGr->cd(5);
2361   fgrIterMeanTheta=new TGraph(1);
2362   fgrIterRMSTheta=new TGraph(1);
2363   fgrIterRMSTheta->GetYaxis()->SetRangeUser(-1000.,1000.);
2364   fgrIterRMSTheta->SetName("fgrIterRMSTheta");
2365   fgrIterRMSTheta->SetLineColor(2);
2366   fgrIterMeanTheta->SetName("fgrIterMeanTheta");
2367   fgrIterMeanTheta->SetTitle("Convergence of #deltaTheta \n");
2368   fgrIterMeanTheta->GetXaxis()->SetTitle("mdeg");
2369   fgrIterRMSTheta->Draw("acp");
2370   fgrIterMeanTheta->Draw("cp");
2371
2372   fCanvGr->cd(6);
2373   fgrIterMeanPhi=new TGraph(1);
2374   fgrIterRMSPhi=new TGraph(1);
2375   fgrIterRMSPhi->GetYaxis()->SetRangeUser(-1000.,1000.);
2376   fgrIterRMSPhi->SetName("fgrIterRMSPhi");
2377   fgrIterRMSPhi->SetLineColor(2);
2378   fgrIterMeanPhi->SetName("fgrIterMeanPhi");
2379   fgrIterMeanPhi->SetTitle("Convergence of #deltaPhi \n");
2380   fgrIterMeanPhi->GetXaxis()->SetTitle("mdeg");
2381   fgrIterRMSPhi->Draw("acp");
2382   fgrIterMeanPhi->Draw("cp");
2383
2384
2385
2386 }
2387
2388 void AliITSRealignTracks::UpdateDraw(TArrayI *volids,Int_t iter,Int_t color){
2389
2390
2391   TString name="hX_";
2392   name+=iter;
2393   name.Append("iter");
2394   TH1F *hX=new TH1F("hX",name.Data(),1000,-10000.,10000.);
2395
2396   name="hY_";
2397   name+=iter;
2398   name.Append("iter");
2399   TH1F *hY=new TH1F("hY",name.Data(),1000,-10000.,10000.);
2400
2401   name="hZ_";
2402   name+=iter;
2403   name.Append("iter");
2404   TH1F *hZ=new TH1F("hZ",name.Data(),1000,-10000.,10000.);
2405
2406   name="hPsi_";
2407   name+=iter;
2408   name.Append("iter");
2409   TH1F *hPsi=new TH1F("hPsi",name.Data(),1000,-5000.,5000.);
2410
2411   name="hTheta_";
2412   name+=iter;
2413   name.Append("iter");
2414   TH1F *hTheta=new TH1F("hTheta",name.Data(),1000,-5000.,5000.);
2415
2416   name="hPhi_";
2417   name+=iter;
2418   name.Append("iter");
2419   TH1F *hPhi=new TH1F("hPhi",name.Data(),1000,-5000.,5000.);
2420
2421   Int_t layer,mod;
2422   Double_t transl[3],rot[3],transldr[3],rotdr[3];
2423
2424   for(Int_t i=0;i<volids->GetSize();i++){
2425     layer=AliGeomManager::VolUIDToLayer(volids->At(i),mod);
2426     fAlignObjs[layer-AliGeomManager::kFirstLayer][mod]->GetPars(transl,rot);
2427     fAlignDrawObjs[layer-AliGeomManager::kFirstLayer][mod]->GetPars(transldr,rotdr);
2428
2429     hX->Fill(10000.*(transl[0]-transldr[0]));
2430     hY->Fill(10000.*(transl[1]-transldr[1]));
2431     hZ->Fill(10000.*(transl[2]-transldr[2]));
2432     hPsi->Fill(1000.*(rot[0]-rotdr[0]));
2433     hTheta->Fill(1000.*(rot[1]-rotdr[1]));
2434     hPhi->Fill(1000.*(rot[1]-rotdr[2]));
2435     //Update the pars of the draw object
2436     fAlignDrawObjs[layer-AliGeomManager::kFirstLayer][mod]->SetPars(transl[0],transl[1],transl[2],rot[0],rot[1],rot[2]);
2437   }
2438
2439   hX->SetLineColor(color);
2440   hY->SetLineColor(color);
2441   hZ->SetLineColor(color);
2442   hPsi->SetLineColor(color);
2443   hTheta->SetLineColor(color);
2444   hPhi->SetLineColor(color);
2445
2446
2447   fCanvPar->cd(1);
2448   hX->Draw("Same");
2449   fCanvPar->cd(2);
2450   hY->Draw("Same");
2451   fCanvPar->cd(3);
2452   hZ->Draw("Same");
2453   fCanvPar->cd(4);
2454   hPsi->Draw("Same");
2455   fCanvPar->cd(5);
2456   hTheta->Draw("Same");
2457   fCanvPar->cd(6);
2458   hPhi->Draw("Same");
2459   gPad->Modified();
2460   fCanvPar->Update();
2461   fCanvPar->Modified();
2462
2463   fgrIterMeanX->SetPoint(fgrIterMeanX->GetN()+1,iter,hX->GetMean());
2464   fgrIterRMSX->SetPoint(fgrIterRMSX->GetN()+1,iter,hX->GetRMS());
2465   fgrIterMeanY->SetPoint(fgrIterMeanY->GetN()+1,iter,hY->GetMean());
2466   fgrIterRMSY->SetPoint(fgrIterRMSY->GetN()+1,iter,hY->GetRMS());
2467   fgrIterMeanZ->SetPoint(fgrIterMeanZ->GetN()+1,iter,hZ->GetMean());
2468   fgrIterRMSZ->SetPoint(fgrIterRMSZ->GetN()+1,iter,hZ->GetRMS());
2469   fgrIterMeanPsi->SetPoint(fgrIterMeanPsi->GetN()+1,iter,hPsi->GetMean());
2470   fgrIterRMSPsi->SetPoint(fgrIterRMSPsi->GetN()+1,iter,hPsi->GetRMS());
2471   fgrIterMeanTheta->SetPoint(fgrIterMeanTheta->GetN()+1,iter,hTheta->GetMean());
2472   fgrIterRMSTheta->SetPoint(fgrIterRMSTheta->GetN()+1,iter,hTheta->GetRMS());
2473   fgrIterMeanPhi->SetPoint(fgrIterMeanPhi->GetN()+1,iter,hPhi->GetMean());
2474   fgrIterRMSPhi->SetPoint(fgrIterRMSPhi->GetN()+1,iter,hPhi->GetRMS());
2475
2476   gPad->Modified();
2477   fCanvGr->Update();
2478   fCanvGr->Update();
2479 }
2480
2481 void AliITSRealignTracks::WriteHists(const char *outfile){
2482
2483   TFile *f=new TFile(outfile,"RECREATE");
2484   f->cd();
2485   fCanvPar->Write();
2486   fCanvGr->Write();
2487   fgrIterMeanX->Write();
2488   fgrIterRMSX->Write();
2489   fgrIterMeanY->Write();
2490   fgrIterRMSY->Write();
2491   fgrIterMeanZ->Write();
2492   fgrIterRMSZ->Write();
2493   fgrIterMeanPsi->Write();
2494   fgrIterRMSPsi->Write();
2495   fgrIterMeanTheta->Write();
2496   fgrIterRMSTheta->Write();
2497   fgrIterMeanPhi->Write();
2498   fgrIterRMSPhi->Write();
2499
2500   f->Close();
2501   return;
2502
2503 }