ITS/AliITSResidualsAnalysis.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 //   Implementation of the tracks residuals analysis class
  18 //   It provides an access to the track space points
  19 //   written along the esd tracks. The class enables
  20 //   the user to plug any track fitter (deriving from
  21 //   AliTrackFitter class) and minimization fo the
  22 //   track residual sums (deriving from the AliTrackResiduals).
  23 //-----------------------------------------------------------------
  24
  25 #include <Riostream.h>
  26 #include <TFile.h>
  27 #include <TMath.h>
  28 #include <TArrayI.h>
  29 #include <TClonesArray.h>
  30 #include <TNtuple.h>
  31 #include <TTree.h>
  32 #include <TF1.h>
  33 #include <TH1F.h>
  34 #include <TH2F.h>
  35 #include <TCanvas.h>
  36 #include <TGraphErrors.h>
  37 #include "TGeoMatrix.h"
  38 #include "TGeoManager.h"
  39 #include "TGeoPhysicalNode.h"
  40 #include "TMatrixDSym.h"
  41 #include "TMatrixDSymEigen.h"
  42 #include "TMatrixD.h"
  43 #include "TString.h"
  44
  45 #include "AliAlignmentTracks.h"
  46 #include "AliTrackPointArray.h"
  47 #include "AliAlignObjParams.h"
  48 #include "AliTrackResiduals.h"
  49 #include "AliTrackFitter.h"
  50 #include "AliTrackFitterKalman.h"
  51 #include "AliTrackFitterRieman.h"
  52 #include "AliTrackResiduals.h"
  53 #include "AliTrackResidualsChi2.h"
  54 #include "AliTrackResidualsFast.h"
  55 #include "AliLog.h"
  56 #include "AliITSgeomTGeo.h"
  57
  58 #include "AliITSResidualsAnalysis.h"
  59
  60 ClassImp(AliITSResidualsAnalysis)
  61
  62 //____________________________________________________________________________
  63   AliITSResidualsAnalysis::AliITSResidualsAnalysis():
  64     AliAlignmentTracks(),
  65     fnHist(0),
  66     fnPhi(0),
  67     fnZ(0),
  68     fvolidsToBin(0),
  69     fLastVolVolid(0),
  70     fCoordToBinTable(0),
  71     fVolResHistRPHI(0),
  72     fResHistZ(0),
  73     fResHistX(0),
  74     fResHistXLocsddL(0),
  75     fResHistXLocsddR(0),
  76     fHistCoordGlobY(0),
  77     fPullHistRPHI(0),
  78     fPullHistZ(0),
  79     fTrackDirPhi(0),
  80     fTrackDirLambda(0),
  81     fTrackDirLambda2(0),
  82     fTrackDirAlpha(0),
  83     fTrackDirPhiAll(0),
  84     fTrackDirLambdaAll(0),
  85     fTrackDirLambda2All(0),
  86     fTrackDirAlphaAll(0),
  87     fTrackDir(0),
  88     fTrackDirAll(0),
  89     fTrackDir2All(0),
  90     fTrackDirXZAll(0),
  91     fResHistGlob(0),
  92     fhistCorrVol(0),
  93     fVolNTracks(0),
  94     fhEmpty(0),
  95     fhistVolNptsUsed(0),
  96     fhistVolUsed(0),
  97     fSigmaVolZ(0),
  98     fsingleLayer(0),
  99     fWriteHist(0),
 100     fpTrackVolIDs(0),
 101     fVolVolids(0),
 102     fVolUsed(0),
 103     fRealignObjFileIsOpen(kFALSE),
 104     fClonesArray(0),
 105     fAliTrackPoints("AliTrackPoints.root"),
 106     fGeom("geometry.root")
 107 {
 108
 109   //
 110   // Defaults
 111   //
 112
 113
 114 }
 115
 116 //____________________________________________________________________________
 117 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TString aliTrackPoints,
 118                                                  const TString geom):
 119   AliAlignmentTracks(),
 120   fnHist(0),
 121   fnPhi(0),
 122   fnZ(0),
 123   fvolidsToBin(0),
 124   fLastVolVolid(0),
 125   fCoordToBinTable(0),
 126   fVolResHistRPHI(0),
 127   fResHistZ(0),
 128   fResHistX(0),
 129   fResHistXLocsddL(0),
 130   fResHistXLocsddR(0),
 131   fHistCoordGlobY(0),
 132   fPullHistRPHI(0),
 133   fPullHistZ(0),
 134   fTrackDirPhi(0),
 135   fTrackDirLambda(0),
 136   fTrackDirLambda2(0),
 137   fTrackDirAlpha(0),
 138   fTrackDirPhiAll(0),
 139   fTrackDirLambdaAll(0),
 140   fTrackDirLambda2All(0),
 141   fTrackDirAlphaAll(0),
 142   fTrackDir(0),
 143   fTrackDirAll(0),
 144   fTrackDir2All(0),
 145   fTrackDirXZAll(0),
 146   fResHistGlob(0),
 147   fhistCorrVol(0),
 148   fVolNTracks(0),
 149   fhEmpty(0),
 150   fhistVolNptsUsed(0),
 151   fhistVolUsed(0),
 152   fSigmaVolZ(0),
 153   fsingleLayer(0),
 154   fWriteHist(0),
 155   fpTrackVolIDs(0),
 156   fVolVolids(0),
 157   fVolUsed(0),
 158   fRealignObjFileIsOpen(kFALSE),
 159   fClonesArray(0),
 160   fAliTrackPoints(aliTrackPoints),
 161   fGeom(geom)
 162 {
 163   //
 164   // Standard Constructor (alitrackpoints)
 165   //
 166
 167
 168 }
 169
 170 //____________________________________________________________________________
 171 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TArrayI *volIDs):
 172   AliAlignmentTracks(),
 173   fnHist(0),
 174   fnPhi(0),
 175   fnZ(0),
 176   fvolidsToBin(0),
 177   fLastVolVolid(0),
 178   fCoordToBinTable(0),
 179   fVolResHistRPHI(0),
 180   fResHistZ(0),
 181   fResHistX(0),
 182   fResHistXLocsddL(0),
 183   fResHistXLocsddR(0),
 184   fHistCoordGlobY(0),
 185   fPullHistRPHI(0),
 186   fPullHistZ(0),
 187   fTrackDirPhi(0),
 188   fTrackDirLambda(0),
 189   fTrackDirLambda2(0),
 190   fTrackDirAlpha(0),
 191   fTrackDirPhiAll(0),
 192   fTrackDirLambdaAll(0),
 193   fTrackDirLambda2All(0),
 194   fTrackDirAlphaAll(0),
 195   fTrackDir(0),
 196   fTrackDirAll(0),
 197   fTrackDir2All(0),
 198   fTrackDirXZAll(0),
 199   fResHistGlob(0),
 200   fhistCorrVol(0),
 201   fVolNTracks(0),
 202   fhEmpty(0),
 203   fhistVolNptsUsed(0),
 204   fhistVolUsed(0),
 205   fSigmaVolZ(0),
 206   fsingleLayer(0),
 207   fWriteHist(0),
 208   fpTrackVolIDs(0),
 209   fVolVolids(0),
 210   fVolUsed(0),
 211   fRealignObjFileIsOpen(kFALSE),
 212   fClonesArray(0),
 213   fAliTrackPoints("AliTrackPoints.root"),
 214   fGeom("geometry.root")
 215
 216 {
 217   //
 218   // Original Constructor
 219   //
 220   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
 221   InitHistograms(volIDs);
 222
 223 }
 224
 225 //____________________________________________________________________________
 226 AliITSResidualsAnalysis::~AliITSResidualsAnalysis()
 227 {
 228   //
 229   //  Destructor
 230   //
 231
 232   if(fvolidsToBin)        delete[] fvolidsToBin;
 233   if(fLastVolVolid)       delete[] fLastVolVolid;
 234   if(fCoordToBinTable)    delete[] fCoordToBinTable;
 235   if(fHistCoordGlobY)     delete[] fHistCoordGlobY;
 236   if(fVolResHistRPHI)     delete fVolResHistRPHI;
 237   if(fResHistZ)           delete fResHistZ;
 238   if(fResHistX){
 239     for(Int_t i=0; i<fnHist; i++) delete fResHistX[i];
 240     delete [] fResHistX;
 241   }
 242   if(fResHistXLocsddL){
 243     for(Int_t i=0; i<fnHist; i++) delete fResHistXLocsddL[i];
 244     delete [] fResHistXLocsddL;
 245   }
 246   if(fResHistXLocsddR){
 247     for(Int_t i=0; i<fnHist; i++) delete fResHistXLocsddR[i];
 248     delete [] fResHistXLocsddR;
 249   }
 250   if(fPullHistRPHI)       delete fPullHistRPHI;
 251   if(fPullHistZ)          delete fPullHistZ;
 252   if(fTrackDirPhi)        delete fTrackDirPhi;
 253   if(fTrackDirLambda)     delete fTrackDirLambda;
 254   if(fTrackDirLambda2)    delete fTrackDirLambda2;
 255   if(fTrackDirAlpha)      delete fTrackDirAlpha;
 256   if(fTrackDirPhiAll)     delete fTrackDirPhiAll;
 257   if(fTrackDirLambdaAll)  delete fTrackDirLambdaAll;
 258   if(fTrackDirLambda2All) delete fTrackDirLambda2All;
 259   if(fTrackDirAlphaAll)   delete fTrackDirAlphaAll;
 260   if(fTrackDir)           delete fTrackDir;
 261   if(fTrackDirAll)        delete fTrackDirAll;
 262   if(fTrackDir2All)       delete fTrackDir2All;
 263   if(fTrackDirXZAll)      delete fTrackDirXZAll;
 264   if(fResHistGlob)        delete fResHistGlob;
 265   if(fhistCorrVol)        delete fhistCorrVol;
 266   if(fVolNTracks)         delete fVolNTracks;
 267   if(fhEmpty)             delete fhEmpty;
 268   if(fhistVolNptsUsed)    delete fhistVolNptsUsed;
 269   if(fhistVolUsed)        delete fhistVolUsed;
 270   if(fSigmaVolZ)          delete fSigmaVolZ;
 271   if(fpTrackVolIDs)       delete fpTrackVolIDs;
 272   if(fVolVolids)          delete fVolVolids;
 273   if(fVolUsed)            delete fVolUsed;
 274   if(fClonesArray)        delete fClonesArray;
 275
 276   SetFileNameTrackPoints("");
 277   SetFileNameGeometry("");
 278
 279 }
 280
 281 //____________________________________________________________________________
 282 void AliITSResidualsAnalysis::InitHistograms(const TArrayI *volIDs)
 283 {
 284   //
 285   // Method that sets and creates the required hisstrograms
 286   // with the correct binning (it dos not fill them)
 287   //
 288
 289
 290   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
 291
 292   TString histnameRPHI="HistRPHI_volID_",aux;
 293   TString histnameZ="HistZ_volID_";
 294   TString histnameX="HistX_volID_";
 295   TString histnameGlob="HistGlob_volID_";
 296   TString histnameCorrVol="HistCorrVol_volID";
 297   TString histnamePullRPHI="HistPullRPHI_volID_";
 298   TString histnamePullZ="HistPullZ_volID_";
 299
 300   TString histnameDirPhi="HistTrackDirPhi_volID_";
 301   TString histnameDirLambda="HistTrackDirLambda_volID_";
 302   TString histnameDirLambda2="HistTrackDirLambda2_volID_";
 303   TString histnameDirAlpha="HistTrackDirAlpha_volID_";
 304   TString histnameDir="HistTrackDir_volID_";
 305
 306   fnHist=volIDs->GetSize();
 307   fVolResHistRPHI=new TH1F*[fnHist];
 308   fResHistGlob=new TH1F*[fnHist];
 309   fResHistZ=new TH1F*[fnHist];
 310   fResHistX=new TH1F*[fnHist];
 311   fResHistXLocsddL=new TH1F*[fnHist];
 312   fResHistXLocsddR=new TH1F*[fnHist];
 313   fHistCoordGlobY=new TH1F*[fnHist];
 314   fPullHistRPHI=new TH1F*[fnHist];
 315   fPullHistZ=new TH1F*[fnHist];
 316   fhistCorrVol=new TH2F*[fnHist];
 317
 318   fTrackDirPhi=new TH1F*[fnHist];
 319   fTrackDirLambda=new TH1F*[fnHist];
 320   fTrackDirLambda2=new TH1F*[fnHist];
 321   fTrackDirAlpha=new TH1F*[fnHist];
 322
 323   fTrackDirPhiAll=new TH1F("fTrackDirPhiAll","fTrackDirPhiAll",100,-180,180);
 324   fTrackDirLambdaAll=new TH1F("fTrackDirLambdaAll","fTrackDirLambdaAll",100,-180,180);
 325   fTrackDirLambda2All=new TH1F("fTrackDirLambda2All","fTrackDirLambda2All",100,0,180);
 326   fTrackDirAlphaAll=new TH1F("fTrackDirAlphaAll","fTrackDirAlphaAll",100,-180,180);
 327
 328   fTrackDirAll=new TH2F("fTrackDirAll","Hist with trakcs directions",100,-180,180,100,-180,180);
 329   fTrackDir2All=new TH2F("fTrackDir2All","Hist with trakcs directions",100,-180,180,100,-180,180);
 330   fTrackDirXZAll=new TH2F("fTrackDirXZAll","Hist with trakcs directions from XZ ",100,-3.,3.,100,-3.,3.);
 331
 332   fTrackDir=new TH2F*[fnHist];
 333
 334   Bool_t binning;
 335   Float_t **binningZPhi;
 336   Float_t *binningz;
 337   Float_t *binningphi;
 338
 339   binningZPhi=CheckSingleLayer(volIDs);
 340   fvolidsToBin=new Int_t*[fnPhi*fnZ];
 341   binningphi=binningZPhi[0];
 342   binningz=binningZPhi[1];
 343   binning=SetBinning(volIDs,binningphi,binningz);
 344
 345   if(binning){ //ONLY FOR A SINGLE LAYER!
 346     fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
 347     fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",fnPhi,binningphi,fnZ,binningz);
 348     fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
 349     fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",fnPhi,binningphi,fnZ,binningz);
 350     fhEmpty=new TH2F("fhEmpty","Hist for getting binning",fnPhi,binningphi,fnZ,binningz);
 351     fVolNTracks->SetXTitle("Volume #phi");
 352     fVolNTracks->SetYTitle("Volume z ");
 353     fhistVolNptsUsed->SetXTitle("Volume #phi");
 354     fhistVolNptsUsed->SetYTitle("Volume z ");
 355     fhistVolUsed->SetXTitle("Volume #phi");
 356     fhistVolUsed->SetYTitle("Volume z ");
 357     fSigmaVolZ->SetXTitle("Volume #phi");
 358     fSigmaVolZ->SetYTitle("Volume z ");
 359   } else{
 360     fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
 361     fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",50,-3.2,3.2,100,-80,80);
 362     fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
 363     fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",50,-3.2,3.2,100,-80,80);
 364     fhEmpty=new TH2F("fhEmpty","Hist for getting binning",50,-3.2,3.2,100,-80,80);
 365     fVolNTracks->SetXTitle("Volume #phi");
 366     fVolNTracks->SetYTitle("Volume z ");
 367     fhistVolNptsUsed->SetXTitle("Volume #phi");
 368     fhistVolNptsUsed->SetYTitle("Volume z ");
 369     fhistVolUsed->SetXTitle("Volume #phi");
 370     fhistVolUsed->SetYTitle("Volume z ");
 371     fSigmaVolZ->SetXTitle("Volume #phi");
 372     fSigmaVolZ->SetYTitle("Volume z ");
 373   }
 374
 375   fpTrackVolIDs=new TArrayI(fnHist);
 376   fVolUsed=new TArrayI*[fnHist];
 377   fVolVolids=new TArrayI*[fnHist];
 378   fLastVolVolid=new Int_t[fnHist];
 379
 380   for (Int_t nhist=0;nhist<fnHist;nhist++){
 381     fpTrackVolIDs->AddAt(volIDs->At(nhist),nhist);
 382     aux=histnameRPHI;
 383     aux+=volIDs->At(nhist);
 384     fVolResHistRPHI[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
 385     fVolResHistRPHI[nhist]->SetName(aux.Data());
 386     fVolResHistRPHI[nhist]->SetTitle(aux.Data());
 387
 388     aux=histnameZ;
 389     aux+=volIDs->At(nhist);
 390     fResHistZ[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
 391     fResHistZ[nhist]->SetName(aux.Data());
 392     fResHistZ[nhist]->SetTitle(aux.Data());
 393
 394     aux=histnameX;
 395     aux+=volIDs->At(nhist);
 396     fResHistX[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
 397     fResHistX[nhist]->SetName(aux.Data());
 398     fResHistX[nhist]->SetTitle(aux.Data());
 399
 400     aux=histnameX;
 401     aux+=volIDs->At(nhist);
 402     aux.Append("LocalSDDLeft");
 403     fResHistXLocsddL[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
 404     fResHistXLocsddL[nhist]->SetName(aux.Data());
 405     fResHistXLocsddL[nhist]->SetTitle(aux.Data());
 406     aux=histnameX;
 407
 408     aux+=volIDs->At(nhist);
 409     aux.Append("LocalSDDRight");
 410     fResHistXLocsddR[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
 411     fResHistXLocsddR[nhist]->SetName(aux.Data());
 412     fResHistXLocsddR[nhist]->SetTitle(aux.Data());
 413
 414     aux="fHistCoordGlobY";
 415     fHistCoordGlobY[nhist]=new TH1F("histname","histname",24000,-30.,30.);
 416     fHistCoordGlobY[nhist]->SetName(aux.Data());
 417     fHistCoordGlobY[nhist]->SetTitle(aux.Data());
 418
 419     aux=histnamePullRPHI;
 420     aux+=volIDs->At(nhist);
 421     fPullHistRPHI[nhist]=new TH1F("histname","histname",100,-7.,7.);
 422     fPullHistRPHI[nhist]->SetName(aux.Data());
 423     fPullHistRPHI[nhist]->SetTitle(aux.Data());
 424
 425     aux=histnamePullZ;
 426     aux+=volIDs->At(nhist);
 427     fPullHistZ[nhist]=new TH1F("histname","histname",100,-7.,7.);
 428     fPullHistZ[nhist]->SetName(aux.Data());
 429     fPullHistZ[nhist]->SetTitle(aux.Data());
 430
 431     aux=histnameDirPhi;
 432     aux+=volIDs->At(nhist);
 433     fTrackDirPhi[nhist]=new TH1F("histname","histname",100,-180,180);
 434     fTrackDirPhi[nhist]->SetName(aux.Data());
 435     fTrackDirPhi[nhist]->SetTitle(aux.Data());
 436
 437     aux=histnameDirLambda;
 438     aux+=volIDs->At(nhist);
 439     fTrackDirLambda[nhist]=new TH1F("histname","histname",100,0,180);
 440     fTrackDirLambda[nhist]->SetName(aux.Data());
 441     fTrackDirLambda[nhist]->SetTitle(aux.Data());
 442
 443     aux=histnameDirLambda2;
 444     aux+=volIDs->At(nhist);
 445     fTrackDirLambda2[nhist]=new TH1F("histname","histname",100,0,180);
 446     fTrackDirLambda2[nhist]->SetName(aux.Data());
 447     fTrackDirLambda2[nhist]->SetTitle(aux.Data());
 448
 449     aux=histnameDirAlpha;
 450     aux+=volIDs->At(nhist);
 451     fTrackDirAlpha[nhist]=new TH1F("histname","histname",100,-180,180);
 452     fTrackDirAlpha[nhist]->SetName(aux.Data());
 453     fTrackDirAlpha[nhist]->SetTitle(aux.Data());
 454
 455     aux=histnameDir;
 456     aux+=volIDs->At(nhist);
 457     fTrackDir[nhist]=new TH2F("histname","histname",100,-90.,90.,100,-180.,180.);
 458     fTrackDir[nhist]->SetName(aux.Data());
 459     fTrackDir[nhist]->SetTitle(aux.Data());
 460
 461     aux=histnameGlob;
 462     aux+=volIDs->At(nhist);
 463     fResHistGlob[nhist]=new TH1F("histname","histname",400,-0.08,0.08);
 464     fResHistGlob[nhist]->SetName(aux.Data());
 465     fResHistGlob[nhist]->SetTitle(aux.Data());
 466
 467     aux=histnameCorrVol;
 468     aux+=volIDs->At(nhist);
 469     fhistCorrVol[nhist]=new TH2F("histname","histname",50,-3.2,3.2,100,-80,80);
 470
 471
 472     fhistCorrVol[nhist]->SetName(aux.Data());
 473     fhistCorrVol[nhist]->SetTitle(aux.Data());
 474     fhistCorrVol[nhist]->SetXTitle("Volume #varphi");
 475     fhistCorrVol[nhist]->SetYTitle("Volume z ");
 476     fVolVolids[nhist]=new TArrayI(100);
 477     fVolUsed[nhist]=new TArrayI(1000);
 478     fLastVolVolid[nhist]=0;
 479
 480   }
 481   fWriteHist=kFALSE;
 482
 483   return;
 484 }
 485
 486 //____________________________________________________________________________
 487 void AliITSResidualsAnalysis::ListVolUsed(TTree *pointsTree,TArrayI ***arrayIndex,Int_t **lastIndex)
 488 {
 489   //
 490   // This is copied from AliAlignmentClass::LoadPoints() method
 491   //
 492   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
 493   Int_t volIDalignable,volIDpoint,iModule;
 494   AliTrackPoint p;
 495   AliTrackPointArray* array = 0;
 496   pointsTree->SetBranchAddress("SP", &array);
 497
 498
 499   for(Int_t ivol=0;ivol<fnHist;ivol++){
 500     Int_t lastused=0;
 501     volIDalignable=fpTrackVolIDs->At(ivol);
 502     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volIDalignable,iModule);
 503
 504     Int_t nArraysId = lastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 505     TArrayI *index = arrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 506     for (Int_t iArrayId = 0;iArrayId < nArraysId; iArrayId++) {
 507
 508       // Get tree entry
 509       Int_t entry = (*index)[iArrayId];
 510
 511       pointsTree->GetEvent(entry);
 512       if (!array) {
 513         AliWarning("Wrong space point array index!");
 514         continue;
 515       }
 516
 517       // Get the space-point array
 518       Int_t modnum,nPoints = array->GetNPoints();
 519
 520       for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
 521
 522
 523         array->GetPoint(p,iPoint);
 524
 525         AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
 526         // check if the layer id is valid
 527         if ((layer < AliGeomManager::kFirstLayer) ||
 528             (layer >= AliGeomManager::kLastLayer)) {
 529           AliError(Form("Layer index is invalid: %d (%d -> %d) !",
 530                         layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
 531           continue;
 532         }
 533
 534
 535         if ((modnum >= AliGeomManager::LayerSize(layer)) ||
 536             (modnum < 0)) {
 537           AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
 538                         layer,modnum,AliGeomManager::LayerSize(layer)));
 539           continue;
 540         }
 541         if (layer > AliGeomManager::kSSD2) continue; // ITS only
 542
 543
 544         volIDpoint=(Int_t)p.GetVolumeID();
 545         if (volIDpoint==volIDalignable) continue;
 546         Int_t size = fVolVolids[ivol]->GetSize();
 547         // If needed allocate new size
 548         if (fLastVolVolid[ivol]>=size){// Warning: fLAST[NHIST] is useless
 549           fVolVolids[ivol]->Set(size + 1000);
 550         }
 551
 552
 553         fVolVolids[ivol]->AddAt(volIDpoint,fLastVolVolid[ivol]);
 554         fLastVolVolid[ivol]++;
 555         Bool_t usedVol=kFALSE;
 556         for(Int_t used=0;used<lastused;used++){
 557           if(fVolUsed[ivol]->At(used)==volIDpoint){
 558             usedVol=kTRUE;
 559             break;
 560           }
 561         }
 562
 563
 564         if (!usedVol){
 565           size = fVolUsed[ivol]->GetSize();
 566           // If needed allocate new size
 567           if (lastused>= size){
 568             fVolUsed[ivol]->Set(size + 1000);
 569           }
 570           fVolUsed[ivol]->AddAt(volIDpoint,lastused);
 571           lastused++;
 572         }
 573
 574         FillVolumeCorrelationHists(ivol,volIDalignable,volIDpoint,usedVol);
 575
 576       }// end loop
 577     }
 578   }
 579   fWriteHist=kTRUE;
 580   return;
 581 }
 582
 583 //____________________________________________________________________________
 584 void AliITSResidualsAnalysis::FillVolumeCorrelationHists(Int_t ivol,Int_t volIDalignable,Int_t volIDpoint,Bool_t usedVol) const
 585 {
 586   //
 587   // Fill the histograms with the correlations between volumes
 588   //
 589
 590
 591   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
 592   Double_t translGlobal[3];
 593   //  Double_t radius;
 594   Double_t phi;
 595   //  const char *symname,*volpath;
 596   /*  TGeoPNEntry *pne;
 597       TGeoPhysicalNode *pn;
 598       TGeoHMatrix *globMatrix;
 599
 600
 601       symname = AliGeomManager::SymName(volIDalignable);
 602       pne = gGeoManager->GetAlignableEntry(symname);
 603       volpath=pne->GetTitle();
 604       pn=gGeoManager->MakePhysicalNode(volpath);
 605       globMatrix=pn->GetMatrix();
 606   */
 607
 608   AliGeomManager::GetOrigTranslation(volIDalignable,translGlobal);
 609   //  radius=TMath::Sqrt(transGlobal[0]*transGlobal[0]+transGlobal[1]*transGlobal[1]);
 610   phi=TMath::ATan2(translGlobal[1],translGlobal[0]);
 611   fhistVolNptsUsed->Fill(phi,translGlobal[2]);
 612   if(!usedVol){
 613     fhistVolUsed->Fill(phi,translGlobal[2]);
 614   }
 615
 616   /*  symname = AliGeomManager::SymName(volIDpoint);
 617       pne = gGeoManager->GetAlignableEntry(symname);
 618       volpath=pne->GetTitle();
 619       pn=gGeoManager->MakePhysicalNode(volpath);
 620       globMatrix=pn->GetMatrix();
 621       transGlobal=globMatrix->GetTranslation();
 622   */
 623   AliGeomManager::GetOrigTranslation(volIDpoint,translGlobal);
 624   //  radius=TMath::Sqrt(transGlobal[0]*transGlobal[0]+transGlobal[1]*transGlobal[1]);
 625   phi=TMath::ATan2(translGlobal[1],translGlobal[0]);
 626
 627   fhistCorrVol[ivol]->Fill(phi,translGlobal[2]);
 628
 629   return;
 630 }
 631
 632 //____________________________________________________________________________
 633 void AliITSResidualsAnalysis::FillResidualsH(AliTrackPointArray *points,
 634                                              AliTrackPointArray *pTrack) const
 635 {
 636   //
 637   // Method that fills the histograms with the residuals
 638   //
 639
 640   Int_t volIDpoint;
 641   Float_t xyz[3],xyz2[3];
 642   Double_t xyzD[3],xyz2D[3];
 643   Double_t loc[3],loc2[3];
 644
 645   Float_t resRPHI,resGlob,resZ,resX;
 646
 647   Double_t pullrphi,sign,phi;
 648   AliTrackPoint p,pTr;
 649
 650   for(Int_t ipoint=0;ipoint<points->GetNPoints();ipoint++){
 651
 652     //pTrack->GetPoint(pTr,ipoint);
 653     points->GetPoint(p,ipoint);
 654     volIDpoint=(Int_t)p.GetVolumeID();
 655     p.GetXYZ(xyz);
 656
 657     pTrack->GetPoint(pTr,ipoint);
 658     pTr.GetXYZ(xyz2);
 659
 660     for(Int_t i=0;i<3;i++){
 661       xyzD[i]=xyz[i];
 662       xyz2D[i]=xyz2[i];
 663     }
 664
 665     phi = TMath::ATan2(xyz[1],xyz[0]);//<-watch out: phi of the pPoints!
 666
 667     resZ=xyz2[2]-xyz[2];
 668     resX=xyz2[0]-xyz[0];
 669
 670     resRPHI=TMath::Sqrt((xyz2[0]-xyz[0])*(xyz2[0]-xyz[0])+(xyz2[1]-xyz[1])*(xyz2[1]-xyz[1]));
 671
 672     sign=TMath::ATan2(xyz2[1],xyz2[0])-TMath::ATan2(xyz[1],xyz[0]);
 673     if(sign!=0.){
 674       sign=sign/TMath::Abs(sign);
 675       resRPHI=resRPHI*sign;
 676
 677     }
 678     else{
 679       pullrphi=0.;
 680       resRPHI=0.;
 681     }
 682
 683     resGlob=TMath::Sqrt((xyz2[0]-xyz[0])*(xyz2[0]-xyz[0])+(xyz2[1]-xyz[1])*(xyz2[1]-xyz[1])+(xyz2[2]-xyz[2])*(xyz2[2]-xyz[2]));
 684
 685     for(Int_t ivolIDs=0;ivolIDs<fpTrackVolIDs->GetSize();ivolIDs++){
 686       if(volIDpoint==fpTrackVolIDs->At(ivolIDs)){
 687
 688         fVolResHistRPHI[ivolIDs]->Fill(resRPHI);
 689         fResHistZ[ivolIDs]->Fill(resZ);
 690         fResHistX[ivolIDs]->Fill(resX);
 691         fHistCoordGlobY[ivolIDs]->Fill(xyz[1]);
 692
 693         Int_t modIndex = -1; // SDD Section
 694         if(AliGeomManager::VolUIDToLayer(volIDpoint)==3) modIndex=volIDpoint-6144+240;
 695         if(AliGeomManager::VolUIDToLayer(volIDpoint)==4) modIndex=volIDpoint-8192+240+84;
 696         if(modIndex>0){
 697           AliITSgeomTGeo::GlobalToLocal(modIndex,xyzD,loc); // error here!?
 698           AliITSgeomTGeo::GlobalToLocal(modIndex,xyz2D,loc2);
 699           Float_t rexloc=loc2[0]-loc[0];
 700           //cout<<"Residual: "<<volIDpoint<<" "<<loc[0]<<" -> "<<rexloc<<endl;
 701           if(loc[0]<0){
 702             fResHistXLocsddR[ivolIDs]->Fill(rexloc);
 703           }else{
 704             fResHistXLocsddL[ivolIDs]->Fill(rexloc);
 705           }
 706         }
 707         fResHistGlob[ivolIDs]->Fill(resGlob);
 708
 709         fTrackDirPhiAll->Fill(phi);
 710         fTrackDirPhi[ivolIDs]->Fill(phi);
 711
 712         if(fsingleLayer){
 713           Int_t binz,binphi;
 714           Float_t globalPhi,globalZ;
 715           if(kTRUE||(fvolidsToBin[ivolIDs][0]!=volIDpoint)){
 716             binphi=GetBinPhiZ((Int_t)volIDpoint,&binz);
 717           }
 718           else{
 719             // This in the case of alignment of one entire layer
 720             // (fnHIst=layersize) may reduce iterations:
 721             // remind of that fsingleLayer->fnHista<layerSize
 722             binphi=fvolidsToBin[ivolIDs][1];
 723             binz=fvolidsToBin[ivolIDs][2];
 724           }
 725           globalPhi=fCoordToBinTable[binphi][binz][0];
 726           globalZ=fCoordToBinTable[binphi][binz][1];
 727
 728           fVolNTracks->Fill(globalPhi,globalZ);
 729         }
 730         else fVolNTracks->Fill(TMath::ATan2(xyz[1],xyz[0]),xyz[2]);
 731       }
 732     }
 733   }
 734 }
 735
 736 //____________________________________________________________________________
 737 Bool_t AliITSResidualsAnalysis::SaveHists(Int_t minNpoints, TString outname) const
 738 {
 739   //
 740   // Saves the histograms into a tree and saves the tree into a file
 741   //
 742
 743
 744   // Output file
 745   TFile *hFile=new TFile(outname.Data(),"RECREATE","File containing the Residuals Tree");
 746
 747   // TTree with the residuals
 748   TTree *analysisTree=new TTree("analysisTree","Tree with the residuals");
 749
 750   // Declares Variables to be stored into the TTree
 751   TF1 *gauss;
 752   Int_t volID,entries,nHistAnalyzed=0;
 753   Double_t meanResRPHI,meanResZ,meanResX,rmsResRPHI,rmsResZ,rmsResX,coordVol[3],x,y,z;
 754   TH1F *histRPHI = new TH1F();
 755   TH1F *histZ = new TH1F();
 756   TH1F *histX = new TH1F();
 757   TH1F *histXLocsddL = new TH1F();
 758   TH1F *histXLocsddR = new TH1F();
 759   TH1F *histCoordGlobY = new TH1F();
 760   // Note: 0 = RPHI, 1 = Z
 761
 762
 763   // Branching the TTree
 764   analysisTree->Branch("volID",&volID,"volID/I");
 765   analysisTree->Branch("x",&x,"x/D");
 766   analysisTree->Branch("y",&y,"y/D");
 767   analysisTree->Branch("z",&z,"z/D");
 768   analysisTree->Branch("meanResRPHI",&meanResRPHI,"meanResRPHI/D");
 769   analysisTree->Branch("meanResZ",&meanResZ,"meanResZ/D");
 770   analysisTree->Branch("meanResX",&meanResX,"meanResX/D");
 771   analysisTree->Branch("rmsResRPHI",&rmsResRPHI,"rmsResRPHI/D");
 772   analysisTree->Branch("rmsResZ",&rmsResZ,"rmsResZ/D");
 773
 774   analysisTree->Branch("histRPHI","TH1F",&histRPHI,128000,0);
 775   analysisTree->Branch("histZ","TH1F",&histZ,128000,0);
 776   analysisTree->Branch("histX","TH1F",&histX,128000,0);
 777   analysisTree->Branch("histXLocsddL","TH1F",&histXLocsddL,128000,0);
 778   analysisTree->Branch("histXLocsddR","TH1F",&histXLocsddR,128000,0);
 779   analysisTree->Branch("histCoordGlobY","TH1F",&histCoordGlobY,128000,0);
 780
 781   Int_t blimps=0;
 782
 783   for(Int_t j=0;j<fnHist;j++){
 784
 785     volID=fpTrackVolIDs->At(j);
 786     AliGeomManager::GetTranslation(volID,coordVol);
 787     x=coordVol[0];
 788     y=coordVol[1];
 789     z=coordVol[2];
 790
 791     entries=(Int_t)(fResHistGlob[j]->GetEntries());
 792     blimps+=entries;
 793
 794     if(entries>=minNpoints){
 795       nHistAnalyzed++;
 796
 797       // Entries
 798       //entries=(Int_t)fVolResHistRPHI[j]->GetEntries();
 799
 800       // Filling the RPHI
 801       histRPHI=fVolResHistRPHI[j];
 802       rmsResRPHI=fVolResHistRPHI[j]->GetRMS();
 803         // Fit (for average)
 804       gauss=new TF1("gauss","gaus",-3*rmsResRPHI,3*rmsResRPHI);
 805       fVolResHistRPHI[j]->Fit("gauss","QRN");
 806       meanResRPHI=gauss->GetParameter(1);
 807
 808       // Filling the Z
 809       histZ=fResHistZ[j];
 810       rmsResZ=fResHistZ[j]->GetRMS();
 811         // Fit (for average)
 812       gauss=new TF1("gauss","gaus",-3*rmsResZ,3*rmsResZ);
 813       fResHistZ[j]->Fit("gauss","QRN");
 814       meanResZ=gauss->GetParameter(1);
 815
 816       // Filling the X
 817       histX=fResHistX[j];
 818       rmsResX=fResHistX[j]->GetRMS();
 819         // Fit (for average)
 820       gauss=new TF1("gauss","gaus",-3*rmsResX,3*rmsResX);
 821       fResHistX[j]->Fit("gauss","QRN");
 822       meanResX=gauss->GetParameter(1);
 823
 824       histXLocsddL=fResHistXLocsddL[j];
 825       histXLocsddR=fResHistXLocsddR[j];
 826       histCoordGlobY=fHistCoordGlobY[j];
 827
 828       analysisTree->Fill();
 829     }else{
 830
 831       // Entries
 832       //entries=(Int_t)fVolResHistRPHI[j]->GetEntries();
 833
 834       // Filling the RPHI
 835       histRPHI=fVolResHistRPHI[j];
 836       rmsResRPHI=-1.0;
 837       meanResRPHI=0.0;
 838
 839       // Filling the Z
 840       histZ=fResHistZ[j];
 841       rmsResZ=-1.0;
 842       meanResZ=0.0;
 843
 844       // Filling the X
 845       histX=fResHistX[j];
 846       rmsResX=-1.0;
 847       meanResX=0.0;
 848       histXLocsddL=fResHistXLocsddL[j];
 849       histXLocsddR=fResHistXLocsddR[j];
 850       histCoordGlobY=fHistCoordGlobY[j];
 851
 852       analysisTree->Fill();
 853
 854     }
 855
 856   }
 857
 858   cout<<"-> Modules Analyzed: "<<nHistAnalyzed<<endl;
 859   cout<<"   With "<<blimps<<" events"<<endl;
 860
 861   if(blimps>0){
 862     hFile->cd();
 863     analysisTree->Write();
 864     fVolNTracks->Write();
 865     fhEmpty->Write();
 866     if(fWriteHist){
 867       //TCanvas *color = new TCanvas("color","fhistVolUsed",800,600);
 868       //fhistVolUsed->DrawCopy("COLZ");
 869       fSigmaVolZ->Write();
 870       fhistVolUsed->Write();
 871       /*      fTrackDirPhiAll->Write();
 872               fTrackDirLambdaAll->Write();
 873               fTrackDirLambda2All->Write();
 874               fTrackDirAlphaAll->Write();
 875               fTrackDirAll->Write();
 876               fTrackDir2All->Write();
 877               fTrackDirXZAll->Write();
 878               hFile->Close();*/
 879       fhistVolNptsUsed->Write();
 880       hFile->mkdir("CorrVol");
 881       hFile->cd("CorrVol");
 882       for(Int_t corr=0;corr<fnHist;corr++)fhistCorrVol[corr]->Write();
 883     }
 884     hFile->cd();
 885     //    fhistVolNptsUsed->Write();
 886     hFile->Close();
 887     return kTRUE;
 888   }else {
 889     delete analysisTree;
 890     delete hFile;
 891     return kFALSE;}
 892 }
 893
 894 //____________________________________________________________________________
 895 void AliITSResidualsAnalysis::DrawHists() const
 896 {
 897   //
 898   // Draws the histograms of the residuals and of the number of tracks
 899   //
 900
 901   TString cname;
 902   for(Int_t canv=0;canv<fnHist;canv++){
 903     cname="canv Residuals";
 904     cname+=canv;
 905     TCanvas *c=new TCanvas(cname.Data(),cname.Data(),700,700);
 906     c->Divide(3,1);
 907     c->cd(1);
 908     fVolResHistRPHI[canv]->Draw();
 909     c->cd(2);
 910     fResHistZ[canv]->Draw();
 911     c->cd(3);
 912     fResHistGlob[canv]->Draw();
 913   }
 914   cname="canv NVolTracks";
 915
 916   TCanvas *c2=new TCanvas(cname.Data(),cname.Data(),700,700);
 917   c2->cd();
 918   fVolNTracks->Draw();
 919
 920   return;
 921 }
 922
 923 //____________________________________________________________________________
 924 Float_t** AliITSResidualsAnalysis::CheckSingleLayer(const TArrayI *volids)
 925 {
 926   //
 927   // Checks if volumes array is a single (ITS) layer or not
 928   //
 929
 930   Float_t **binningzphi=new Float_t*[2];
 931   Int_t iModule;
 932   AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volids->At(0),iModule);
 933   //Check that one single Layer is going to be aligned
 934   for(Int_t nvol=0;nvol<volids->GetSize();nvol++){
 935     if(iLayer != AliGeomManager::VolUIDToLayer((UShort_t)volids->At(nvol),iModule)){
 936       printf("Wrong Layer! \n %d , %d , %d ,%d \n",(Int_t)AliGeomManager::VolUIDToLayer((UShort_t)volids->At(nvol),iModule),nvol,volids->GetSize(),iModule);
 937       fsingleLayer=kFALSE;
 938       return binningzphi;
 939     }
 940   }
 941
 942   //Bool_t used=kFALSE;
 943   switch (iLayer) {
 944   case AliGeomManager::kSPD1:{
 945     fnPhi=kPhiSPD1;//kPhiSPD1;
 946     fnZ=kZSPD1;//nZSPD1;
 947     binningzphi[0]=new Float_t[kPhiSPD1+1];
 948     binningzphi[1]=new Float_t[kZSPD1+1];
 949     fCoordToBinTable=new Double_t**[kPhiSPD1];
 950     for(Int_t j=0;j<fnPhi;j++){
 951       fCoordToBinTable[j]=new Double_t*[kZSPD1];
 952     }
 953   }; break;
 954   case AliGeomManager::kSPD2:{
 955     fnPhi=kPhiSPD2;//kPhiSPD1;
 956     fnZ=kZSPD2;//nZSPD1;
 957     binningzphi[0]=new Float_t[kPhiSPD2+1];
 958     binningzphi[1]=new Float_t[kZSPD2+1];
 959     fCoordToBinTable=new Double_t**[kPhiSPD2];
 960     for(Int_t j=0;j<fnPhi;j++){
 961       fCoordToBinTable[j]=new Double_t*[kZSPD2];
 962     }
 963
 964   }; break; case AliGeomManager::kSDD1:{
 965     fnPhi=kPhiSDD1;//kPhiSPD1;
 966     fnZ=kZSDD1;//nZSPD1;
 967     binningzphi[0]=new Float_t[kPhiSDD1+1];
 968     binningzphi[1]=new Float_t[kZSDD1+1];
 969     fCoordToBinTable=new Double_t**[kPhiSDD1];
 970     for(Int_t j=0;j<fnPhi;j++){
 971       fCoordToBinTable[j]=new Double_t*[kZSDD1];
 972     }
 973   }; break; case AliGeomManager::kSDD2:{
 974     fnPhi=kPhiSDD2;//kPhiSPD1;
 975     fnZ=kZSDD2;//nZSPD1;
 976     binningzphi[0]=new Float_t[kPhiSDD2+1];
 977     binningzphi[1]=new Float_t[kZSDD2+1];
 978     fCoordToBinTable=new Double_t**[kPhiSDD2];
 979     for(Int_t j=0;j<fnPhi;j++){
 980       fCoordToBinTable[j]=new Double_t*[kZSDD2];
 981     }
 982   }; break; case AliGeomManager::kSSD1:{
 983     fnPhi=kPhiSSD1;//kPhiSPD1;
 984     fnZ=kZSSD1;//nZSPD1;
 985     binningzphi[0]=new Float_t[kPhiSSD1+1];
 986     binningzphi[1]=new Float_t[kZSSD1+1];
 987     fCoordToBinTable=new Double_t**[kPhiSSD1];
 988     for(Int_t j=0;j<fnPhi;j++){
 989       fCoordToBinTable[j]=new Double_t*[kZSSD1];
 990     }
 991   }; break; case AliGeomManager::kSSD2:{
 992     fnPhi=kPhiSSD2;//kPhiSPD1;
 993     fnZ=kZSSD2;//nZSPD1;
 994     binningzphi[0]=new Float_t[kPhiSSD2+1];
 995     binningzphi[1]=new Float_t[kZSSD2+1];
 996     fCoordToBinTable=new Double_t**[kPhiSSD2];
 997     for(Int_t j=0;j<fnPhi;j++){
 998       fCoordToBinTable[j]=new Double_t*[kZSSD2];
 999     }
1000   }; break;
1001
1002   default:{
1003     printf("Wrong Layer Label! \n");
1004     fsingleLayer=kFALSE;
1005     //////////////////////
1006     fnPhi=1;//kPhiSPD1;
1007     fnZ=1;//nZSPD1;
1008     binningzphi[0]=new Float_t[1];
1009     binningzphi[1]=new Float_t[1];
1010     fCoordToBinTable=new Double_t**[1];
1011     for(Int_t j=0;j<fnPhi;j++){
1012       fCoordToBinTable[j]=new Double_t*[1];
1013     }
1014     return binningzphi;
1015     /////////////////////
1016     // return 0x0;
1017   }
1018   }
1019   fsingleLayer=kTRUE;
1020   return binningzphi;
1021 }
1022
1023 //____________________________________________________________________________
1024 Bool_t AliITSResidualsAnalysis::SetBinning(const TArrayI *volids,Float_t *phiBin,Float_t *zBin)
1025 {
1026   //
1027   // Sets the correct binning
1028   //
1029
1030   if(!fsingleLayer)return kFALSE;
1031   //  const char *volpath,*symname;
1032   Int_t iModule;
1033   Int_t *orderArrayPhi,*orderArrayZ;
1034   UShort_t volID;
1035   Double_t *phiArray,*zArray,*phiArrayOrdered,*zArrayOrdered;
1036   Double_t translGlobal[3];
1037   Double_t lastPhimin=-10;
1038   Double_t lastZmin=-99;
1039   Int_t ***orderPhiZ;
1040   /*  TGeoPNEntry *pne;
1041       TGeoPhysicalNode *pn;
1042       TGeoHMatrix *globMatrix;*/
1043
1044   Bool_t used=kFALSE;
1045
1046   orderPhiZ=new Int_t**[fnPhi];
1047   phiArray=new Double_t[fnPhi];//phiBin[nModulesPhi+1];
1048   zArray=new Double_t[fnZ];//zBin[nModulesZ+1];
1049   phiArrayOrdered=new Double_t[fnPhi];
1050   zArrayOrdered=new Double_t[fnZ];
1051   orderArrayPhi=new Int_t[fnPhi];
1052   orderArrayZ=new Int_t[fnZ];
1053   for(Int_t k=0;k<fnZ;k++){
1054     orderArrayZ[k]=0;
1055     zArray[k]=-1000;
1056   }
1057   for(Int_t k=0;k<fnPhi;k++){
1058     orderArrayPhi[k]=0;
1059     phiArray[k]=-10;
1060     orderPhiZ[k]=new Int_t*[fnZ];
1061   }
1062
1063
1064   AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volids->At(0),iModule);
1065
1066   Int_t lastPhi=0,lastZ=0;
1067   for(iModule=0;iModule<AliGeomManager::LayerSize(iLayer);iModule++){
1068     fvolidsToBin[iModule]=new Int_t[3];
1069     volID=AliGeomManager::LayerToVolUID(iLayer,iModule);
1070     fvolidsToBin[iModule][0]=volID;
1071     /*    symname = AliGeomManager::SymName(volID);
1072           pne = gGeoManager->GetAlignableEntry(symname);
1073           volpath=pne->GetTitle();
1074           pn=gGeoManager->MakePhysicalNode(volpath);
1075           globMatrix=pn->GetMatrix();
1076           translGlobal=globMatrix->GetTranslation();
1077
1078     */
1079     AliGeomManager::GetOrigTranslation(volID,translGlobal);
1080
1081     for(Int_t j=0;j<lastPhi;j++){
1082       used=kFALSE;
1083       if(TMath::Abs(phiArray[j]-TMath::ATan2(translGlobal[1],translGlobal[0]))<2*TMath::Pi()/(10*fnPhi)){//10 is a safety factor but....
1084         fvolidsToBin[iModule][1]=j;
1085         used=kTRUE;
1086         break;
1087       }
1088     }
1089     if(!used){
1090       phiArray[lastPhi]=TMath::ATan2(translGlobal[1],translGlobal[0]);
1091       fvolidsToBin[iModule][1]=lastPhi;
1092       if(phiArray[lastPhi]<lastPhimin)lastPhimin=phiArray[lastPhi];
1093       lastPhi++;
1094       if(lastPhi>fnPhi){
1095         printf("Wrong Phi! \n");
1096         return kFALSE;}
1097     }
1098
1099     for(Int_t j=0;j<lastZ;j++){
1100       used=kFALSE;
1101       if(TMath::Abs(zArray[j]-translGlobal[2])<0.1){
1102         fvolidsToBin[iModule][2]=j;
1103         used=kTRUE;
1104         break;
1105       }
1106     }
1107     if(!used){
1108       fvolidsToBin[iModule][2]=lastZ;
1109       zArray[lastZ]=translGlobal[2];
1110       if(zArray[lastZ]<lastZmin)lastZmin=zArray[lastZ];
1111       lastZ++;
1112       if(lastZ>fnZ){
1113         printf("Wrong Z! \n");
1114         return kFALSE;}
1115     }
1116   }
1117
1118
1119   //ORDERING THE ARRAY OF PHI AND Z VALUES
1120   for(Int_t order=0;order<fnPhi;order++){
1121     for(Int_t j=0;j<fnPhi;j++){
1122       if((j!=order)&&(phiArray[order]>phiArray[j]))orderArrayPhi[order]++;
1123     }
1124   }
1125
1126   for(Int_t order=0;order<fnPhi;order++){
1127     for(Int_t j=0;j<fnPhi;j++){
1128       if(orderArrayPhi[j]==order){
1129         phiArrayOrdered[order]=phiArray[j];
1130         break;
1131         }
1132     }
1133   }
1134
1135
1136   for(Int_t order=0;order<fnZ;order++){
1137     for(Int_t j=0;j<fnZ;j++){
1138       if((j!=order)&&(zArray[order]>zArray[j]))orderArrayZ[order]++;
1139     }
1140   }
1141
1142
1143   for(Int_t order=0;order<fnZ;order++){
1144     for(Int_t j=0;j<fnZ;j++){
1145       if(orderArrayZ[j]==order){
1146         zArrayOrdered[order]=zArray[j];
1147         break;
1148       }
1149     }
1150   }
1151
1152
1153   //Filling the  fCoordToBinTable
1154   for(Int_t j=0;j<fnPhi;j++){
1155     for(Int_t i=0;i<fnZ;i++){
1156       orderPhiZ[j][i]=new Int_t[2];
1157       orderPhiZ[j][i][0]=orderArrayPhi[j];
1158       orderPhiZ[j][i][1]=orderArrayZ[i];
1159     }
1160   }
1161
1162
1163   for(Int_t j=0;j<fnPhi;j++){
1164     for(Int_t i=0;i<fnZ;i++){
1165       fCoordToBinTable[j][i]=new Double_t[2];
1166       fCoordToBinTable[j][i][0]=phiArrayOrdered[j];
1167       fCoordToBinTable[j][i][1]=zArrayOrdered[i];
1168       printf("Now (binphi,binz)= %d, %d e (phi,z)=%f,%f \n",j,i,fCoordToBinTable[j][i][0],fCoordToBinTable[j][i][1]);
1169     }
1170   }
1171   Int_t istar,jstar;
1172   for(iModule=0;iModule<fnPhi*fnZ;iModule++){
1173     istar=fvolidsToBin[iModule][1];
1174     jstar=fvolidsToBin[iModule][2];
1175     fvolidsToBin[iModule][1]=orderPhiZ[istar][jstar][0];
1176     fvolidsToBin[iModule][2]=orderPhiZ[istar][jstar][1];
1177   }
1178
1179
1180   //now constructing the binning
1181   for(Int_t iModPhi=0;iModPhi<fnPhi-1;iModPhi++){
1182     phiBin[iModPhi+1]=(Float_t)phiArrayOrdered[iModPhi]+0.5*(phiArrayOrdered[iModPhi+1]-phiArrayOrdered[iModPhi]);
1183   }
1184
1185   phiBin[0]=(Float_t)phiArrayOrdered[0]-(phiArrayOrdered[1]-phiArrayOrdered[0])/2.;
1186
1187   phiBin[fnPhi]=(Float_t)phiArrayOrdered[fnPhi-1]+(phiArrayOrdered[fnPhi-1]-phiArrayOrdered[fnPhi-2])/2.;
1188   for(Int_t iModPhi=0;iModPhi<fnPhi+1;iModPhi++){
1189     printf("Mean Phi mod %d su %d:  %f \n",iModPhi+1,fnPhi,phiBin[iModPhi]);
1190   }
1191
1192   for(Int_t iModZ=0;iModZ<fnZ-1;iModZ++){
1193     zBin[iModZ+1]=(Float_t)zArrayOrdered[iModZ]+0.5*(zArrayOrdered[iModZ+1]-zArrayOrdered[iModZ]);
1194   }
1195   zBin[0]=(Float_t)zArrayOrdered[0]-0.5*(zArrayOrdered[1]-zArrayOrdered[0]);
1196   zBin[fnZ]=(Float_t)zArrayOrdered[fnZ-1]+0.5*(zArrayOrdered[1]-zArrayOrdered[0]);
1197
1198
1199   for(Int_t iModPhi=0;iModPhi<fnZ+1;iModPhi++){
1200      printf("Mean Z mod %d su %d:  %f \n",iModPhi+1,fnZ,zBin[iModPhi]);
1201   }
1202   return kTRUE;
1203 }
1204
1205
1206 //____________________________________________________________________________
1207 Int_t AliITSResidualsAnalysis::GetBinPhiZ(const Int_t volID,Int_t *binz) const
1208 {
1209   //
1210   // Returns the correct Phi-Z bin
1211   //
1212
1213   if (!fsingleLayer){
1214     printf("No Single Layer reAlignment! \n");
1215     return 100;
1216   }
1217
1218   for(Int_t j=0;j<fnPhi*fnZ;j++){
1219     if(j==fnZ*fnPhi){
1220       printf("Wrong volume UID introduced! fnHist: %d  volID: %d iter: %d \n",fnHist,volID,j);
1221       return 100;
1222     }
1223     if(fvolidsToBin[j][0]==volID){
1224
1225       *binz=fvolidsToBin[j][2];
1226       return fvolidsToBin[j][1];
1227     }
1228   }
1229
1230   return 100;
1231
1232 }
1233
1234 //___________________________________________________________________________
1235 Int_t AliITSResidualsAnalysis::WhichSector(Int_t module) const
1236 {
1237   //
1238   // This method returns the number of the SPD Sector
1239   // to which belongs the module (Sectors 0-9)
1240
1241   //--->cSect = 0 <---
1242   if(module==2048
1243      || module==2049
1244      || module==2050
1245      || module==2051
1246      || module==2052
1247      || module==2053
1248      || module==2054
1249      || module==2055
1250      || module==4096
1251      || module==4097
1252      || module==4098
1253      || module==4099
1254      || module==4100
1255      || module==4101
1256      || module==4102
1257      || module==4103
1258      || module==4104
1259      || module==4105
1260      || module==4106
1261      || module==4107
1262      || module==4108
1263      || module==4109
1264      || module==4110
1265      || module==4111) return 0;
1266
1267   //--->cSect = 1 <---
1268   if(module==2056
1269      || module==2057
1270      || module==2058
1271      || module==2059
1272      || module==2060
1273      || module==2061
1274      || module==2062
1275      || module==2063
1276      || module==4112
1277      || module==4113
1278      || module==4114
1279      || module==4115
1280      || module==4116
1281      || module==4117
1282      || module==4118
1283      || module==4119
1284      || module==4120
1285      || module==4121
1286      || module==4122
1287      || module==4123
1288      || module==4124
1289      || module==4125
1290      || module==4126
1291      || module==4127) return 1;
1292
1293   //--->cSect = 2 <---
1294   if(module==2064
1295      || module==2065
1296      || module==2066
1297      || module==2067
1298      || module==2068
1299      || module==2069
1300      || module==2070
1301      || module==2071
1302      || module==4128
1303      || module==4129
1304      || module==4130
1305      || module==4131
1306      || module==4132
1307      || module==4133
1308      || module==4134
1309      || module==4135
1310      || module==4136
1311      || module==4137
1312      || module==4138
1313      || module==4139
1314      || module==4140
1315      || module==4141
1316      || module==4142
1317      || module==4143) return 2;
1318
1319   //--->cSect = 3 <---
1320   if(module==2072
1321      || module==2073
1322      || module==2074
1323      || module==2075
1324      || module==2076
1325      || module==2077
1326      || module==2078
1327      || module==2079
1328      || module==4144
1329      || module==4145
1330      || module==4146
1331      || module==4147
1332      || module==4148
1333      || module==4149
1334      || module==4150
1335      || module==4151
1336      || module==4152
1337      || module==4153
1338      || module==4154
1339      || module==4155
1340      || module==4156
1341      || module==4157
1342      || module==4158
1343      || module==4159) return 3;
1344
1345   //--->cSect = 4 <---
1346   if(module==2080
1347      || module==2081
1348      || module==2082
1349      || module==2083
1350      || module==2084
1351      || module==2085
1352      || module==2086
1353      || module==2087
1354      || module==4160
1355      || module==4161
1356      || module==4162
1357      || module==4163
1358      || module==4164
1359      || module==4165
1360      || module==4166
1361      || module==4167
1362      || module==4168
1363      || module==4169
1364      || module==4170
1365      || module==4171
1366      || module==4172
1367      || module==4173
1368      || module==4174
1369      || module==4175) return 4;
1370
1371   //--->cSect = 5 <---
1372   if(module==2088
1373      || module==2089
1374      || module==2090
1375      || module==2091
1376      || module==2092
1377      || module==2093
1378      || module==2094
1379      || module==2095
1380      || module==4176
1381      || module==4177
1382      || module==4178
1383      || module==4179
1384      || module==4180
1385      || module==4181
1386      || module==4182
1387      || module==4183
1388      || module==4184
1389      || module==4185
1390      || module==4186
1391      || module==4187
1392      || module==4188
1393      || module==4189
1394      || module==4190
1395      || module==4191) return 5;
1396
1397   //--->cSect = 6 <---
1398   if(module==2096
1399      || module==2097
1400      || module==2098
1401      || module==2099
1402      || module==2100
1403      || module==2101
1404      || module==2102
1405      || module==2103
1406      || module==4192
1407      || module==4193
1408      || module==4194
1409      || module==4195
1410      || module==4196
1411      || module==4197
1412      || module==4198
1413      || module==4199
1414      || module==4200
1415      || module==4201
1416      || module==4202
1417      || module==4203
1418      || module==4204
1419      || module==4205
1420      || module==4206
1421      || module==4207) return 6;
1422
1423   //--->cSect = 7 <---
1424   if(module==2104
1425      || module==2105
1426      || module==2106
1427      || module==2107
1428      || module==2108
1429      || module==2109
1430      || module==2110
1431      || module==2111
1432      || module==4208
1433      || module==4209
1434      || module==4210
1435      || module==4211
1436      || module==4212
1437      || module==4213
1438      || module==4214
1439      || module==4215
1440      || module==4216
1441      || module==4217
1442      || module==4218
1443      || module==4219
1444      || module==4220
1445      || module==4221
1446      || module==4222
1447      || module==4223) return 7;
1448
1449   //--->cSect = 8 <---
1450   if(module==2112
1451      || module==2113
1452      || module==2114
1453      || module==2115
1454      || module==2116
1455      || module==2117
1456      || module==2118
1457      || module==2119
1458      || module==4224
1459      || module==4225
1460      || module==4226
1461      || module==4227
1462      || module==4228
1463      || module==4229
1464      || module==4230
1465      || module==4231
1466      || module==4232
1467      || module==4233
1468      || module==4234
1469      || module==4235
1470      || module==4236
1471      || module==4237
1472      || module==4238
1473      || module==4239) return 8;
1474
1475   //--->cSect = 9 <---
1476   if(module==2120
1477      || module==2121
1478      || module==2122
1479      || module==2123
1480      || module==2124
1481      || module==2125
1482      || module==2126
1483      || module==2127
1484      || module==4240
1485      || module==4241
1486      || module==4242
1487      || module==4243
1488      || module==4244
1489      || module==4245
1490      || module==4246
1491      || module==4247
1492      || module==4248
1493      || module==4249
1494      || module==4250
1495      || module==4251
1496      || module==4252
1497      || module==4253
1498      || module==4254
1499      || module==4255) return 9;
1500
1501   //printf("Module not belonging to SPD, sorry!");
1502   return -1;
1503
1504 }
1505
1506 //____________________________________________________________________________
1507 TArrayI* AliITSResidualsAnalysis::GetSPDSectorsVolids(Int_t sectors[10]) const
1508 {
1509   //
1510   // This method gets the volID Array for the chosen sectors.
1511   // You have to pass an array with a 1 for each selected sector.
1512   // i.e. sectors[10] = {1,1,0,0,0,0,0,0,1,0} -> Sector 0, 1, 9 selected.
1513   //
1514
1515   Int_t nSect=0;
1516   Int_t iModule=0;
1517
1518   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1519
1520   for(Int_t co=0;co<10;co++){ //counts the number of sectors chosen
1521     if(sectors[co]==1) nSect++;
1522   }
1523
1524   if(nSect<1){ //if no sector chosen -> exit
1525     Printf("Error! No Sector/s Selected!");
1526     return 0x0;
1527   }
1528
1529   TArrayI *volIDs = new TArrayI(nSect*24);
1530
1531     if(sectors[0]==1){ //--->cSect = 0 <---
1532       volIDs->AddAt(2048,iModule); iModule++;
1533       volIDs->AddAt(2049,iModule); iModule++;
1534       volIDs->AddAt(2050,iModule); iModule++;
1535       volIDs->AddAt(2051,iModule); iModule++;
1536       volIDs->AddAt(2052,iModule); iModule++;
1537       volIDs->AddAt(2053,iModule); iModule++;
1538       volIDs->AddAt(2054,iModule); iModule++;
1539       volIDs->AddAt(2055,iModule); iModule++;
1540       volIDs->AddAt(4096,iModule); iModule++;
1541       volIDs->AddAt(4097,iModule); iModule++;
1542       volIDs->AddAt(4098,iModule); iModule++;
1543       volIDs->AddAt(4099,iModule); iModule++;
1544       volIDs->AddAt(4100,iModule); iModule++;
1545       volIDs->AddAt(4101,iModule); iModule++;
1546       volIDs->AddAt(4102,iModule); iModule++;
1547       volIDs->AddAt(4103,iModule); iModule++;
1548       volIDs->AddAt(4104,iModule); iModule++;
1549       volIDs->AddAt(4105,iModule); iModule++;
1550       volIDs->AddAt(4106,iModule); iModule++;
1551       volIDs->AddAt(4107,iModule); iModule++;
1552       volIDs->AddAt(4108,iModule); iModule++;
1553       volIDs->AddAt(4109,iModule); iModule++;
1554       volIDs->AddAt(4110,iModule); iModule++;
1555       volIDs->AddAt(4111,iModule); iModule++;
1556     }
1557     if(sectors[1]==1){ //--->cSect = 1 <//---
1558       volIDs->AddAt(2056,iModule); iModule++;
1559       volIDs->AddAt(2057,iModule); iModule++;
1560       volIDs->AddAt(2058,iModule); iModule++;
1561       volIDs->AddAt(2059,iModule); iModule++;
1562       volIDs->AddAt(2060,iModule); iModule++;
1563       volIDs->AddAt(2061,iModule); iModule++;
1564       volIDs->AddAt(2062,iModule); iModule++;
1565       volIDs->AddAt(2063,iModule); iModule++;
1566       volIDs->AddAt(4112,iModule); iModule++;
1567       volIDs->AddAt(4113,iModule); iModule++;
1568       volIDs->AddAt(4114,iModule); iModule++;
1569       volIDs->AddAt(4115,iModule); iModule++;
1570       volIDs->AddAt(4116,iModule); iModule++;
1571       volIDs->AddAt(4117,iModule); iModule++;
1572       volIDs->AddAt(4118,iModule); iModule++;
1573       volIDs->AddAt(4119,iModule); iModule++;
1574       volIDs->AddAt(4120,iModule); iModule++;
1575       volIDs->AddAt(4121,iModule); iModule++;
1576       volIDs->AddAt(4122,iModule); iModule++;
1577       volIDs->AddAt(4123,iModule); iModule++;
1578       volIDs->AddAt(4124,iModule); iModule++;
1579       volIDs->AddAt(4125,iModule); iModule++;
1580       volIDs->AddAt(4126,iModule); iModule++;
1581       volIDs->AddAt(4127,iModule); iModule++;
1582     }
1583     if(sectors[2]==1){//--->cSect = 2 <//---
1584       volIDs->AddAt(2064,iModule); iModule++;
1585       volIDs->AddAt(2065,iModule); iModule++;
1586       volIDs->AddAt(2066,iModule); iModule++;
1587       volIDs->AddAt(2067,iModule); iModule++;
1588       volIDs->AddAt(2068,iModule); iModule++;
1589       volIDs->AddAt(2069,iModule); iModule++;
1590       volIDs->AddAt(2070,iModule); iModule++;
1591       volIDs->AddAt(2071,iModule); iModule++;
1592       volIDs->AddAt(4128,iModule); iModule++;
1593       volIDs->AddAt(4129,iModule); iModule++;
1594       volIDs->AddAt(4130,iModule); iModule++;
1595       volIDs->AddAt(4131,iModule); iModule++;
1596       volIDs->AddAt(4132,iModule); iModule++;
1597       volIDs->AddAt(4133,iModule); iModule++;
1598       volIDs->AddAt(4134,iModule); iModule++;
1599       volIDs->AddAt(4135,iModule); iModule++;
1600       volIDs->AddAt(4136,iModule); iModule++;
1601       volIDs->AddAt(4137,iModule); iModule++;
1602       volIDs->AddAt(4138,iModule); iModule++;
1603       volIDs->AddAt(4139,iModule); iModule++;
1604       volIDs->AddAt(4140,iModule); iModule++;
1605       volIDs->AddAt(4141,iModule); iModule++;
1606       volIDs->AddAt(4142,iModule); iModule++;
1607       volIDs->AddAt(4143,iModule); iModule++;
1608     }
1609     if(sectors[3]==1){//--->cSect = 3 <//---
1610       volIDs->AddAt(2072,iModule); iModule++;
1611       volIDs->AddAt(2073,iModule); iModule++;
1612       volIDs->AddAt(2074,iModule); iModule++;
1613       volIDs->AddAt(2075,iModule); iModule++;
1614       volIDs->AddAt(2076,iModule); iModule++;
1615       volIDs->AddAt(2077,iModule); iModule++;
1616       volIDs->AddAt(2078,iModule); iModule++;
1617       volIDs->AddAt(2079,iModule); iModule++;
1618       volIDs->AddAt(4144,iModule); iModule++;
1619       volIDs->AddAt(4145,iModule); iModule++;
1620       volIDs->AddAt(4146,iModule); iModule++;
1621       volIDs->AddAt(4147,iModule); iModule++;
1622       volIDs->AddAt(4148,iModule); iModule++;
1623       volIDs->AddAt(4149,iModule); iModule++;
1624       volIDs->AddAt(4150,iModule); iModule++;
1625       volIDs->AddAt(4151,iModule); iModule++;
1626       volIDs->AddAt(4152,iModule); iModule++;
1627       volIDs->AddAt(4153,iModule); iModule++;
1628       volIDs->AddAt(4154,iModule); iModule++;
1629       volIDs->AddAt(4155,iModule); iModule++;
1630       volIDs->AddAt(4156,iModule); iModule++;
1631       volIDs->AddAt(4157,iModule); iModule++;
1632       volIDs->AddAt(4158,iModule); iModule++;
1633       volIDs->AddAt(4159,iModule); iModule++;
1634     }
1635     if(sectors[4]==1){//--->cSect = 4 <//---
1636       volIDs->AddAt(2080,iModule); iModule++;
1637       volIDs->AddAt(2081,iModule); iModule++;
1638       volIDs->AddAt(2082,iModule); iModule++;
1639       volIDs->AddAt(2083,iModule); iModule++;
1640       volIDs->AddAt(2084,iModule); iModule++;
1641       volIDs->AddAt(2085,iModule); iModule++;
1642       volIDs->AddAt(2086,iModule); iModule++;
1643       volIDs->AddAt(2087,iModule); iModule++;
1644       volIDs->AddAt(4160,iModule); iModule++;
1645       volIDs->AddAt(4161,iModule); iModule++;
1646       volIDs->AddAt(4162,iModule); iModule++;
1647       volIDs->AddAt(4163,iModule); iModule++;
1648       volIDs->AddAt(4164,iModule); iModule++;
1649       volIDs->AddAt(4165,iModule); iModule++;
1650       volIDs->AddAt(4166,iModule); iModule++;
1651       volIDs->AddAt(4167,iModule); iModule++;
1652       volIDs->AddAt(4168,iModule); iModule++;
1653       volIDs->AddAt(4169,iModule); iModule++;
1654       volIDs->AddAt(4170,iModule); iModule++;
1655       volIDs->AddAt(4171,iModule); iModule++;
1656       volIDs->AddAt(4172,iModule); iModule++;
1657       volIDs->AddAt(4173,iModule); iModule++;
1658       volIDs->AddAt(4174,iModule); iModule++;
1659       volIDs->AddAt(4175,iModule); iModule++;
1660     }
1661     if(sectors[5]==1){//--->cSect = 5 <//---
1662       volIDs->AddAt(2088,iModule); iModule++;
1663       volIDs->AddAt(2089,iModule); iModule++;
1664       volIDs->AddAt(2090,iModule); iModule++;
1665       volIDs->AddAt(2091,iModule); iModule++;
1666       volIDs->AddAt(2092,iModule); iModule++;
1667       volIDs->AddAt(2093,iModule); iModule++;
1668       volIDs->AddAt(2094,iModule); iModule++;
1669       volIDs->AddAt(2095,iModule); iModule++;
1670       volIDs->AddAt(4176,iModule); iModule++;
1671       volIDs->AddAt(4177,iModule); iModule++;
1672       volIDs->AddAt(4178,iModule); iModule++;
1673       volIDs->AddAt(4179,iModule); iModule++;
1674       volIDs->AddAt(4180,iModule); iModule++;
1675       volIDs->AddAt(4181,iModule); iModule++;
1676       volIDs->AddAt(4182,iModule); iModule++;
1677       volIDs->AddAt(4183,iModule); iModule++;
1678       volIDs->AddAt(4184,iModule); iModule++;
1679       volIDs->AddAt(4185,iModule); iModule++;
1680       volIDs->AddAt(4186,iModule); iModule++;
1681       volIDs->AddAt(4187,iModule); iModule++;
1682       volIDs->AddAt(4188,iModule); iModule++;
1683       volIDs->AddAt(4189,iModule); iModule++;
1684       volIDs->AddAt(4190,iModule); iModule++;
1685       volIDs->AddAt(4191,iModule); iModule++;
1686     }
1687     if(sectors[6]==1){//--->cSect = 6 <//---
1688       volIDs->AddAt(2096,iModule); iModule++;
1689       volIDs->AddAt(2097,iModule); iModule++;
1690       volIDs->AddAt(2098,iModule); iModule++;
1691       volIDs->AddAt(2099,iModule); iModule++;
1692       volIDs->AddAt(2100,iModule); iModule++;
1693       volIDs->AddAt(2101,iModule); iModule++;
1694       volIDs->AddAt(2102,iModule); iModule++;
1695       volIDs->AddAt(2103,iModule); iModule++;
1696       volIDs->AddAt(4192,iModule); iModule++;
1697       volIDs->AddAt(4193,iModule); iModule++;
1698       volIDs->AddAt(4194,iModule); iModule++;
1699       volIDs->AddAt(4195,iModule); iModule++;
1700       volIDs->AddAt(4196,iModule); iModule++;
1701       volIDs->AddAt(4197,iModule); iModule++;
1702       volIDs->AddAt(4198,iModule); iModule++;
1703       volIDs->AddAt(4199,iModule); iModule++;
1704       volIDs->AddAt(4200,iModule); iModule++;
1705       volIDs->AddAt(4201,iModule); iModule++;
1706       volIDs->AddAt(4202,iModule); iModule++;
1707       volIDs->AddAt(4203,iModule); iModule++;
1708       volIDs->AddAt(4204,iModule); iModule++;
1709       volIDs->AddAt(4205,iModule); iModule++;
1710       volIDs->AddAt(4206,iModule); iModule++;
1711       volIDs->AddAt(4207,iModule); iModule++;
1712     }
1713      if(sectors[7]==1){ //--->cSect = 7 <//---
1714        volIDs->AddAt(2104,iModule); iModule++;
1715        volIDs->AddAt(2105,iModule); iModule++;
1716        volIDs->AddAt(2106,iModule); iModule++;
1717        volIDs->AddAt(2107,iModule); iModule++;
1718        volIDs->AddAt(2108,iModule); iModule++;
1719        volIDs->AddAt(2109,iModule); iModule++;
1720        volIDs->AddAt(2110,iModule); iModule++;
1721        volIDs->AddAt(2111,iModule); iModule++;
1722        volIDs->AddAt(4208,iModule); iModule++;
1723        volIDs->AddAt(4209,iModule); iModule++;
1724        volIDs->AddAt(4210,iModule); iModule++;
1725        volIDs->AddAt(4211,iModule); iModule++;
1726        volIDs->AddAt(4212,iModule); iModule++;
1727        volIDs->AddAt(4213,iModule); iModule++;
1728        volIDs->AddAt(4214,iModule); iModule++;
1729        volIDs->AddAt(4215,iModule); iModule++;
1730        volIDs->AddAt(4216,iModule); iModule++;
1731        volIDs->AddAt(4217,iModule); iModule++;
1732        volIDs->AddAt(4218,iModule); iModule++;
1733        volIDs->AddAt(4219,iModule); iModule++;
1734        volIDs->AddAt(4220,iModule); iModule++;
1735        volIDs->AddAt(4221,iModule); iModule++;
1736        volIDs->AddAt(4222,iModule); iModule++;
1737        volIDs->AddAt(4223,iModule); iModule++;
1738      }
1739      if(sectors[8]==1){//--->cSect = 8 <//---
1740        volIDs->AddAt(2112,iModule); iModule++;
1741        volIDs->AddAt(2113,iModule); iModule++;
1742        volIDs->AddAt(2114,iModule); iModule++;
1743        volIDs->AddAt(2115,iModule); iModule++;
1744        volIDs->AddAt(2116,iModule); iModule++;
1745        volIDs->AddAt(2117,iModule); iModule++;
1746        volIDs->AddAt(2118,iModule); iModule++;
1747        volIDs->AddAt(2119,iModule); iModule++;
1748        volIDs->AddAt(4224,iModule); iModule++;
1749        volIDs->AddAt(4225,iModule); iModule++;
1750        volIDs->AddAt(4226,iModule); iModule++;
1751        volIDs->AddAt(4227,iModule); iModule++;
1752        volIDs->AddAt(4228,iModule); iModule++;
1753        volIDs->AddAt(4229,iModule); iModule++;
1754        volIDs->AddAt(4230,iModule); iModule++;
1755        volIDs->AddAt(4231,iModule); iModule++;
1756        volIDs->AddAt(4232,iModule); iModule++;
1757        volIDs->AddAt(4233,iModule); iModule++;
1758        volIDs->AddAt(4234,iModule); iModule++;
1759        volIDs->AddAt(4235,iModule); iModule++;
1760        volIDs->AddAt(4236,iModule); iModule++;
1761        volIDs->AddAt(4237,iModule); iModule++;
1762        volIDs->AddAt(4238,iModule); iModule++;
1763        volIDs->AddAt(4239,iModule); iModule++;
1764      }
1765      if(sectors[9]==1){//--->cSect = 9 <//---
1766        volIDs->AddAt(2120,iModule); iModule++;
1767        volIDs->AddAt(2121,iModule); iModule++;
1768        volIDs->AddAt(2122,iModule); iModule++;
1769        volIDs->AddAt(2123,iModule); iModule++;
1770        volIDs->AddAt(2124,iModule); iModule++;
1771        volIDs->AddAt(2125,iModule); iModule++;
1772        volIDs->AddAt(2126,iModule); iModule++;
1773        volIDs->AddAt(2127,iModule); iModule++;
1774        volIDs->AddAt(4240,iModule); iModule++;
1775        volIDs->AddAt(4241,iModule); iModule++;
1776        volIDs->AddAt(4242,iModule); iModule++;
1777        volIDs->AddAt(4243,iModule); iModule++;
1778        volIDs->AddAt(4244,iModule); iModule++;
1779        volIDs->AddAt(4245,iModule); iModule++;
1780        volIDs->AddAt(4246,iModule); iModule++;
1781        volIDs->AddAt(4247,iModule); iModule++;
1782        volIDs->AddAt(4248,iModule); iModule++;
1783        volIDs->AddAt(4249,iModule); iModule++;
1784        volIDs->AddAt(4250,iModule); iModule++;
1785        volIDs->AddAt(4251,iModule); iModule++;
1786        volIDs->AddAt(4252,iModule); iModule++;
1787        volIDs->AddAt(4253,iModule); iModule++;
1788        volIDs->AddAt(4254,iModule); iModule++;
1789        volIDs->AddAt(4255,iModule); iModule++;
1790      }
1791
1792   return volIDs;
1793
1794 }
1795
1796 //____________________________________________________________________________
1797 TArrayI* AliITSResidualsAnalysis::GetITSLayersVolids(Int_t layers[6]) const
1798 {
1799   //
1800   // This method gets the volID Array for the chosen layers.
1801   // You have to pass an array with a 1 for each selected layer.
1802   // i.e. layers[6] = {1,1,0,0,1,1} -> SPD + SSD
1803   //
1804
1805   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1806
1807   Int_t size=0,last=0;
1808
1809   // evaluates the size of the array
1810   for(Int_t i=0;i<6;i++) if(layers[i]==1) size+=AliGeomManager::LayerSize(i+1);
1811
1812   if(size==0){
1813     printf("Error: no layer selected");
1814     return 0x0;
1815   }
1816
1817   TArrayI *volids = new TArrayI(size);
1818
1819   // fills the volId array only for the chosen layers
1820   for(Int_t ilayer=1;ilayer<7;ilayer++){
1821
1822     if(layers[ilayer-1]!=1) continue;
1823
1824     for(Int_t imod=0;imod<AliGeomManager::LayerSize(ilayer);imod++){
1825       volids->AddAt(AliGeomManager::LayerToVolUID(ilayer,imod),last);
1826       last++;
1827     }
1828   }
1829
1830   return volids;
1831
1832 }
1833
1834 //____________________________________________________________________________
1835 void AliITSResidualsAnalysis::GetTrackDirClusterCov(AliTrackPoint *point,Double_t &phi,Double_t &lambda,Double_t &lambda2,Double_t &alpha,Double_t &xovery,Double_t &zovery) const
1836 {
1837   //
1838   // ...
1839   //
1840
1841   TMatrixDSym cov(3);
1842   const Float_t *covvector=point->GetCov();
1843   cov(0,0)=covvector[0];
1844   cov(1,0)=cov(0,1)=covvector[1];
1845   cov(2,0)=cov(0,2)=covvector[2];
1846   cov(1,1)=covvector[3];
1847   cov(1,2)=cov(2,1)=covvector[4];
1848   cov(2,2)=covvector[5];
1849
1850   Double_t determinant=cov.Determinant();
1851   if(determinant!=0.){
1852     TMatrixD vect(3,3);
1853     TVectorD eigenvalues(3);
1854     const TMatrixDSymEigen keigen(cov);
1855     eigenvalues=keigen.GetEigenValues();
1856     vect=keigen.GetEigenVectors();
1857     Double_t mainvect[3];
1858     mainvect[0]=vect(0,0);
1859     mainvect[1]=vect(1,0);
1860     mainvect[2]=vect(2,0);
1861     if(mainvect[1]!=0.){
1862       xovery=mainvect[0]/mainvect[1];
1863       zovery=mainvect[2]/mainvect[1];
1864     }
1865     else {
1866       xovery=9999.;
1867       zovery=9999.;
1868     }
1869     if(mainvect[1]<0.){
1870       mainvect[0]=-1.*mainvect[0];
1871       mainvect[1]=-1.*mainvect[1];
1872       mainvect[2]=-1.*mainvect[2];
1873     }
1874     lambda2=TMath::ATan2(TMath::Sqrt(mainvect[0]*mainvect[0]+mainvect[2]*mainvect[2]),mainvect[1])*TMath::RadToDeg();
1875     lambda=TMath::ATan2(mainvect[2],TMath::Sqrt(mainvect[0]*mainvect[0]+mainvect[1]*mainvect[1]))*TMath::RadToDeg();
1876     phi=TMath::ATan2(mainvect[0],mainvect[2])*TMath::RadToDeg();
1877     alpha=TMath::ATan2(mainvect[1],mainvect[0])*TMath::RadToDeg();
1878   }
1879   else printf("determinant =0!, skip this point \n");
1880
1881   return;
1882 }
1883
1884 //____________________________________________________________________________
1885 void AliITSResidualsAnalysis::CalculateResiduals(const TArrayI *volids,
1886                                                  const TArrayI *volidsfit,
1887                                                  AliGeomManager::ELayerID layerRangeMin,
1888                                                  AliGeomManager::ELayerID layerRangeMax,
1889                                                  TString outname)
1890 {
1891   // CalculateResiduals for a set of detector volumes.
1892   // Tracks are fitted only within
1893   // the range defined by the user
1894   // (by layerRangeMin and layerRangeMax)
1895   // or within the set of volidsfit
1896   // Repeat the procedure 'iterations' times
1897
1898   Int_t nVolIds = volids->GetSize();
1899   if (nVolIds == 0) { AliError("Volume IDs array is empty!"); return; }
1900
1901   // Load only the tracks with at least one
1902   // space point in the set of volume (volids)
1903
1904
1905   //AliAlignmentTracks::SetPointsFilename(GetFileNameTrackPoints());
1906   AliAlignmentTracks::BuildIndex();
1907
1908   ListVolUsed(fPointsTree,fArrayIndex,fLastIndex);
1909   AliTrackPointArray **points;
1910
1911   Int_t pointsDim;
1912   LoadPoints(volids, points,pointsDim);
1913
1914   Int_t nArrays = fPointsTree->GetEntries();
1915
1916   if (nArrays == 0){ AliError("Points array is empty!"); return; }
1917   AliTrackFitter *fitter = CreateFitter();
1918
1919   Int_t ecount=0;
1920   Int_t totcount=0;
1921
1922   for (Int_t iArray = 0; iArray < nArrays; iArray++){
1923
1924     cout<<"Investigating "<<iArray<<"/"<<nArrays<<endl;
1925
1926     if (!points[iArray]){
1927       cout<<" Skipping: "<<iArray<<endl;
1928       continue;
1929     }
1930
1931     fitter->SetTrackPointArray(points[iArray],kTRUE); // Watch out, problems
1932                                                       // when few sectors
1933
1934     totcount++;
1935
1936     // *** FITTING ***
1937     if(fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE){
1938       ecount++;
1939       cout<<"->BAD: "<<iArray<<endl;
1940       continue;
1941     } //else cout<<"->GOOD: "<<iArray<<endl;
1942
1943     AliTrackPointArray *pVolId,*pTrack;
1944
1945     fitter->GetTrackResiduals(pVolId,pTrack);
1946     FillResidualsH(pVolId,pTrack);
1947
1948   }
1949
1950   cout<<"   -> nVolIds: "<<nVolIds<<endl;
1951   cout<<"   -> Non-Fitted tracks: "<<ecount<<"/"<<totcount<<endl;
1952
1953   UnloadPoints(pointsDim,points);
1954
1955   SaveHists(3,outname);
1956
1957   return;
1958
1959 }
1960
1961
1962 //______________________________________________________________________________
1963 void AliITSResidualsAnalysis::ProcessVolumes(Int_t fit,
1964                                              TArrayI *volIDs,
1965                                              TArrayI *volIDsFit,
1966                                              TString misalignmentFile,
1967                                              TString outname,
1968                                              Int_t minPoints)
1969 {
1970   //
1971   // This function process the AliTrackPoints and volID (into residuals)
1972   //
1973
1974   // setting up geometry and the trackpoints file
1975   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1976
1977   SetPointsFilename(GetFileNameTrackPoints());
1978
1979   // creating some tools
1980   AliTrackFitter *fitter;
1981   if(fit==1){
1982     fitter = new AliTrackFitterKalman();
1983   }else fitter = new AliTrackFitterRieman();
1984
1985   fitter->SetMinNPoints(minPoints);
1986
1987   SetTrackFitter(fitter);
1988
1989   if(misalignmentFile=="")printf("NO FAKE MISALIGNMENT\n");
1990   else {
1991     Bool_t misal=Misalign(misalignmentFile,"ITSAlignObjs");
1992     if(!misal){
1993       printf("PROBLEM WITH FAKE MISALIGNMENT!");
1994       return;
1995     }
1996   }
1997
1998   CalculateResiduals(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,outname);
1999
2000     return;
2001
2002 }