ITS/AliITSResidualsAnalysis.cxx

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