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