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