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 <TChain.h>
  27 #include <TFile.h>
  28 #include <TMath.h>
  29 #include <TArrayI.h>
  30 #include <TClonesArray.h>
  31 #include <TNtuple.h>
  32 #include <TTree.h>
  33 #include <TF1.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TCanvas.h>
  37 #include <TGraphErrors.h>
  38 #include "TGeoMatrix.h"
  39 #include "TGeoManager.h"
  40 #include "TGeoPhysicalNode.h"
  41 #include "TMatrixDSym.h"
  42 #include "TMatrixDSymEigen.h"
  43 #include "TMatrixD.h"
  44 #include "TString.h"
  45
  46 #include "AliAlignmentTracks.h"
  47 #include "AliTrackPointArray.h"
  48 #include "AliAlignObjParams.h"
  49 #include "AliTrackResiduals.h"
  50 #include "AliTrackFitter.h"
  51 #include "AliTrackFitterKalman.h"
  52 #include "AliTrackFitterRieman.h"
  53 #include "AliTrackResiduals.h"
  54 #include "AliTrackResidualsChi2.h"
  55 #include "AliTrackResidualsFast.h"
  56 #include "AliLog.h"
  57
  58 #include "AliITSResidualsAnalysis.h"
  59
  60
  61 // Structure of the RealignmentAnalysisLayer.C
  62   typedef struct {
  63     Int_t nentries;
  64     Float_t rms;
  65     Float_t meanFit;
  66     Float_t errmeanFit;
  67     Float_t sigmaFit;
  68   }  histProperties_t;
  69
  70 ClassImp(AliITSResidualsAnalysis)
  71
  72 //____________________________________________________________________________
  73   AliITSResidualsAnalysis::AliITSResidualsAnalysis():
  74     AliAlignmentTracks(),
  75     fnHist(0),
  76     fnPhi(0),
  77     fnZ(0),
  78     fvolidsToBin(0),
  79     fLastVolVolid(0),
  80     fCoordToBinTable(0),
  81     fVolResHistRPHI(0),
  82     fResHistZ(0),
  83     fPullHistRPHI(0),
  84     fPullHistZ(0),
  85     fTrackDirPhi(0),
  86     fTrackDirLambda(0),
  87     fTrackDirLambda2(0),
  88     fTrackDirAlpha(0),
  89     fTrackDirPhiAll(0),
  90     fTrackDirLambdaAll(0),
  91     fTrackDirLambda2All(0),
  92     fTrackDirAlphaAll(0),
  93     fTrackDir(0),
  94     fTrackDirAll(0),
  95     fTrackDir2All(0),
  96     fTrackDirXZAll(0),
  97     fResHistGlob(0),
  98     fhistCorrVol(0),
  99     fVolNTracks(0),
 100     fhEmpty(0),
 101     fhistVolNptsUsed(0),
 102     fhistVolUsed(0),
 103     fSigmaVolZ(0),
 104     fsingleLayer(0),
 105     fWriteHist(0),
 106     fpTrackVolIDs(0),
 107     fVolVolids(0),
 108     fVolUsed(0),
 109     fRealignObjFileIsOpen(kFALSE),
 110     fClonesArray(0),
 111     fAliTrackPoints("AliTrackPoints.root"),
 112     fGeom("geometry.root")
 113
 114 {
 115
 116   //
 117   // Defaults
 118   //
 119
 120
 121 }
 122
 123 //____________________________________________________________________________
 124 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TString aliTrackPoints,
 125                                                  const TString geom):
 126   AliAlignmentTracks(),
 127   fnHist(0),
 128   fnPhi(0),
 129   fnZ(0),
 130   fvolidsToBin(0),
 131   fLastVolVolid(0),
 132   fCoordToBinTable(0),
 133   fVolResHistRPHI(0),
 134   fResHistZ(0),
 135   fPullHistRPHI(0),
 136   fPullHistZ(0),
 137   fTrackDirPhi(0),
 138   fTrackDirLambda(0),
 139   fTrackDirLambda2(0),
 140   fTrackDirAlpha(0),
 141   fTrackDirPhiAll(0),
 142   fTrackDirLambdaAll(0),
 143   fTrackDirLambda2All(0),
 144   fTrackDirAlphaAll(0),
 145   fTrackDir(0),
 146   fTrackDirAll(0),
 147   fTrackDir2All(0),
 148   fTrackDirXZAll(0),
 149   fResHistGlob(0),
 150   fhistCorrVol(0),
 151   fVolNTracks(0),
 152   fhEmpty(0),
 153   fhistVolNptsUsed(0),
 154   fhistVolUsed(0),
 155   fSigmaVolZ(0),
 156   fsingleLayer(0),
 157   fWriteHist(0),
 158   fpTrackVolIDs(0),
 159   fVolVolids(0),
 160   fVolUsed(0),
 161   fRealignObjFileIsOpen(kFALSE),
 162   fClonesArray(0),
 163   fAliTrackPoints(aliTrackPoints),
 164   fGeom(geom)
 165 {
 166   //
 167   // Constructor (alitrackpoints)
 168   //
 169
 170
 171 }
 172
 173 //____________________________________________________________________________
 174 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TArrayI *volIDs):
 175   AliAlignmentTracks(),
 176   fnHist(0),
 177   fnPhi(0),
 178   fnZ(0),
 179   fvolidsToBin(0),
 180   fLastVolVolid(0),
 181   fCoordToBinTable(0),
 182   fVolResHistRPHI(0),
 183   fResHistZ(0),
 184   fPullHistRPHI(0),
 185   fPullHistZ(0),
 186   fTrackDirPhi(0),
 187   fTrackDirLambda(0),
 188   fTrackDirLambda2(0),
 189   fTrackDirAlpha(0),
 190   fTrackDirPhiAll(0),
 191   fTrackDirLambdaAll(0),
 192   fTrackDirLambda2All(0),
 193   fTrackDirAlphaAll(0),
 194   fTrackDir(0),
 195   fTrackDirAll(0),
 196   fTrackDir2All(0),
 197   fTrackDirXZAll(0),
 198   fResHistGlob(0),
 199   fhistCorrVol(0),
 200   fVolNTracks(0),
 201   fhEmpty(0),
 202   fhistVolNptsUsed(0),
 203   fhistVolUsed(0),
 204   fSigmaVolZ(0),
 205   fsingleLayer(0),
 206   fWriteHist(0),
 207   fpTrackVolIDs(0),
 208   fVolVolids(0),
 209   fVolUsed(0),
 210   fRealignObjFileIsOpen(kFALSE),
 211   fClonesArray(0),
 212   fAliTrackPoints("AliTrackPoints.root"),
 213   fGeom("geometry.root")
 214
 215 {
 216   //
 217   // Original Constructor
 218   //
 219
 220   InitHistograms(volIDs);
 221
 222 }
 223
 224 //____________________________________________________________________________
 225 AliITSResidualsAnalysis::AliITSResidualsAnalysis(TArrayI *volIDs,AliTrackPointArray **tracksClustArray,AliTrackPointArray **tracksFitPointsArray):
 226   AliAlignmentTracks(),
 227   fnHist(0),
 228   fnPhi(0),
 229   fnZ(0),
 230   fvolidsToBin(0),
 231   fLastVolVolid(0),
 232   fCoordToBinTable(0),
 233   fVolResHistRPHI(0),
 234   fResHistZ(0),
 235   fPullHistRPHI(0),
 236   fPullHistZ(0),
 237   fTrackDirPhi(0),
 238   fTrackDirLambda(0),
 239   fTrackDirLambda2(0),
 240   fTrackDirAlpha(0),
 241   fTrackDirPhiAll(0),
 242   fTrackDirLambdaAll(0),
 243   fTrackDirLambda2All(0),
 244   fTrackDirAlphaAll(0),
 245   fTrackDir(0),
 246   fTrackDirAll(0),
 247   fTrackDir2All(0),
 248   fTrackDirXZAll(0),
 249   fResHistGlob(0),
 250   fhistCorrVol(0),
 251   fVolNTracks(0),
 252   fhEmpty(0),
 253   fhistVolNptsUsed(0),
 254   fhistVolUsed(0),
 255   fSigmaVolZ(0),
 256   fsingleLayer(0),
 257   fWriteHist(0),
 258   fpTrackVolIDs(0),
 259   fVolVolids(0),
 260   fVolUsed(0),
 261   fRealignObjFileIsOpen(kFALSE),
 262   fClonesArray(0),
 263   fAliTrackPoints("AliTrackPoints.root"),
 264   fGeom("geometry.root")
 265
 266 {
 267   //
 268   // Detailed Constructor (deprecated)
 269   //
 270
 271
 272   TString histnameRPHI="HistRPHI_volID_",aux;
 273   TString histnameZ="HistZ_volID_";
 274   TString histnameGlob="HistGlob_volID_";
 275   TString histnameCorrVol="HistCorrVol_volID";
 276   TString histnamePullRPHI="HistPullRPHI_volID_";
 277   TString histnamePullZ="HistPullZ_volID_";
 278   fnHist=volIDs->GetSize();
 279   fVolResHistRPHI=new TH1F*[fnHist];
 280   fResHistGlob=new TH1F*[fnHist];
 281   fResHistZ=new TH1F*[fnHist];
 282   fPullHistRPHI=new TH1F*[fnHist];
 283   fPullHistZ=new TH1F*[fnHist];
 284   fhistCorrVol=new TH2F*[fnHist];
 285   Float_t **binningZPhi=CheckSingleLayer(volIDs);
 286   fvolidsToBin=new Int_t*[fnPhi*fnZ];
 287   Float_t *binningphi=binningZPhi[0];
 288   Float_t *binningz=binningZPhi[1];
 289   Bool_t binning=SetBinning(volIDs,binningphi,binningz);
 290   if(binning){
 291     fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
 292     fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",fnPhi,binningphi,fnZ,binningz);
 293     fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
 294     fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",fnPhi,binningphi,fnZ,binningz);
 295     fhEmpty=new TH2F("fhEmpty","Hist for getting binning",fnPhi,binningphi,fnZ,binningz);
 296     fVolNTracks->SetXTitle("Volume #phi");
 297     fVolNTracks->SetYTitle("Volume z ");
 298     fhistVolNptsUsed->SetXTitle("Volume #phi");
 299     fhistVolNptsUsed->SetYTitle("Volume z ");
 300     fhistVolUsed->SetXTitle("Volume #phi");
 301     fhistVolUsed->SetYTitle("Volume z ");
 302     fSigmaVolZ->SetXTitle("Volume #phi");
 303     fSigmaVolZ->SetYTitle("Volume z ");
 304   }
 305   else{
 306     fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
 307     fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",50,-3.2,3.2,100,-80,80);
 308     fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
 309     fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",50,-3.2,3.2,100,-80,80);
 310     fhEmpty=new TH2F("fhEmpty","Hist for getting binning",50,-3.2,3.2,100,-80,80);
 311     fVolNTracks->SetXTitle("Volume #phi");
 312     fVolNTracks->SetYTitle("Volume z ");
 313     fhistVolNptsUsed->SetXTitle("Volume #phi");
 314     fhistVolNptsUsed->SetYTitle("Volume z ");
 315     fhistVolUsed->SetXTitle("Volume #phi");
 316     fhistVolUsed->SetYTitle("Volume z ");
 317     fSigmaVolZ->SetXTitle("Volume #phi");
 318     fSigmaVolZ->SetYTitle("Volume z ");
 319   }
 320
 321   fpTrackVolIDs=new TArrayI(fnHist);
 322   fVolUsed=new TArrayI*[fnHist];
 323   fVolVolids=new TArrayI*[fnHist];
 324   fLastVolVolid=new Int_t[fnHist];
 325
 326   for (Int_t nhist=0;nhist<fnHist;nhist++){
 327     fpTrackVolIDs->AddAt(volIDs->At(nhist),nhist);
 328     aux=histnameRPHI;
 329     aux+=volIDs->At(nhist);
 330     fVolResHistRPHI[nhist]=new TH1F("namehist","histname",200,-0.02,0.02);
 331     fVolResHistRPHI[nhist]->SetName(aux.Data());
 332     fVolResHistRPHI[nhist]->SetTitle(aux.Data());
 333
 334     aux=histnameZ;
 335     aux+=volIDs->At(nhist);
 336     fResHistZ[nhist]=new TH1F("histname","histname",400,-0.08,0.08);
 337     fResHistZ[nhist]->SetName(aux.Data());
 338     fResHistZ[nhist]->SetTitle(aux.Data());
 339
 340     aux=histnamePullRPHI;
 341     aux+=volIDs->At(nhist);
 342     fPullHistRPHI[nhist]=new TH1F("histname","histname",100,-7.,7.);
 343     fPullHistRPHI[nhist]->SetName(aux.Data());
 344     fPullHistRPHI[nhist]->SetTitle(aux.Data());
 345
 346     aux=histnamePullZ;
 347     aux+=volIDs->At(nhist);
 348     fPullHistZ[nhist]=new TH1F("histname","histname",100,-7.,7.);
 349     fPullHistZ[nhist]->SetName(aux.Data());
 350     fPullHistZ[nhist]->SetTitle(aux.Data());
 351
 352     aux=histnameGlob;
 353     aux+=volIDs->At(nhist);
 354     fResHistGlob[nhist]=new TH1F("histname","histname",400,-0.08,0.08);
 355     fResHistGlob[nhist]->SetName(aux.Data());
 356     fResHistGlob[nhist]->SetTitle(aux.Data());
 357
 358     aux=histnameCorrVol;
 359     aux+=volIDs->At(nhist);
 360     fhistCorrVol[nhist]=new TH2F("histname","histname",50,-3.2,3.2,100,-80,80);
 361     fhistCorrVol[nhist]->SetName(aux.Data());
 362     fhistCorrVol[nhist]->SetTitle(aux.Data());
 363     fhistCorrVol[nhist]->SetXTitle("Volume #varphi");
 364     fhistCorrVol[nhist]->SetYTitle("Volume z ");
 365
 366     fVolVolids[nhist]=new TArrayI(1000);
 367     fVolUsed[nhist]=new TArrayI(1000);
 368     fLastVolVolid[nhist]=0;
 369     FillResHists(tracksClustArray[nhist],tracksFitPointsArray[nhist]);
 370   }
 371   fWriteHist=kFALSE;
 372   DrawHists();
 373
 374   SetFileNameTrackPoints("");  // Filename with the AliTrackPoints
 375   SetFileNameGeometry(""); // Filename with the Geometry
 376
 377
 378 }
 379
 380 //____________________________________________________________________________
 381 AliITSResidualsAnalysis::~AliITSResidualsAnalysis()
 382 {
 383   //
 384   //  Destructor
 385   //
 386
 387   if(fvolidsToBin)        delete[] fvolidsToBin;
 388   if(fLastVolVolid)       delete[] fLastVolVolid;
 389   if(fCoordToBinTable)    delete[] fCoordToBinTable;
 390   if(fVolResHistRPHI)     delete fVolResHistRPHI;
 391   if(fResHistZ)           delete fResHistZ;
 392   if(fPullHistRPHI)       delete fPullHistRPHI;
 393   if(fPullHistZ)          delete fPullHistZ;
 394   if(fTrackDirPhi)        delete fTrackDirPhi;
 395   if(fTrackDirLambda)     delete fTrackDirLambda;
 396   if(fTrackDirLambda2)    delete fTrackDirLambda2;
 397   if(fTrackDirAlpha)      delete fTrackDirAlpha;
 398   if(fTrackDirPhiAll)     delete fTrackDirPhiAll;
 399   if(fTrackDirLambdaAll)  delete fTrackDirLambdaAll;
 400   if(fTrackDirLambda2All) delete fTrackDirLambda2All;
 401   if(fTrackDirAlphaAll)   delete fTrackDirAlphaAll;
 402   if(fTrackDir)           delete fTrackDir;
 403   if(fTrackDirAll)        delete fTrackDirAll;
 404   if(fTrackDir2All)       delete fTrackDir2All;
 405   if(fTrackDirXZAll)      delete fTrackDirXZAll;
 406   if(fResHistGlob)        delete fResHistGlob;
 407   if(fhistCorrVol)        delete fhistCorrVol;
 408   if(fVolNTracks)         delete fVolNTracks;
 409   if(fhEmpty)             delete fhEmpty;
 410   if(fhistVolNptsUsed)    delete fhistVolNptsUsed;
 411   if(fhistVolUsed)        delete fhistVolUsed;
 412   if(fSigmaVolZ)          delete fSigmaVolZ;
 413   if(fpTrackVolIDs)       delete fpTrackVolIDs;
 414   if(fVolVolids)          delete fVolVolids;
 415   if(fVolUsed)            delete fVolUsed;
 416   if(fClonesArray)        delete fClonesArray;
 417
 418   SetFileNameTrackPoints("");
 419   SetFileNameGeometry("");
 420
 421 }
 422
 423 //____________________________________________________________________________
 424 void AliITSResidualsAnalysis::InitHistograms(const TArrayI *volIDs)
 425 {
 426   //
 427   // Method that sets and creates the required hisstrograms
 428   // with the correct binning (it dos not fill them)
 429   //
 430
 431   TString histnameRPHI="HistRPHI_volID_",aux;
 432   TString histnameZ="HistZ_volID_";
 433   TString histnameGlob="HistGlob_volID_";
 434   TString histnameCorrVol="HistCorrVol_volID";
 435   TString histnamePullRPHI="HistPullRPHI_volID_";
 436   TString histnamePullZ="HistPullZ_volID_";
 437
 438   TString histnameDirPhi="HistTrackDirPhi_volID_";
 439   TString histnameDirLambda="HistTrackDirLambda_volID_";
 440   TString histnameDirLambda2="HistTrackDirLambda2_volID_";
 441   TString histnameDirAlpha="HistTrackDirAlpha_volID_";
 442   TString histnameDir="HistTrackDir_volID_";
 443
 444
 445   fnHist=volIDs->GetSize();
 446   fVolResHistRPHI=new TH1F*[fnHist];
 447   fResHistGlob=new TH1F*[fnHist];
 448   fResHistZ=new TH1F*[fnHist];
 449   fPullHistRPHI=new TH1F*[fnHist];
 450   fPullHistZ=new TH1F*[fnHist];
 451   fhistCorrVol=new TH2F*[fnHist];
 452
 453
 454   fTrackDirPhi=new TH1F*[fnHist];
 455   fTrackDirLambda=new TH1F*[fnHist];
 456   fTrackDirLambda2=new TH1F*[fnHist];
 457   fTrackDirAlpha=new TH1F*[fnHist];
 458
 459
 460   fTrackDirPhiAll=new TH1F("fTrackDirPhiAll","fTrackDirPhiAll",100,-180,180);
 461   fTrackDirLambdaAll=new TH1F("fTrackDirLambdaAll","fTrackDirLambdaAll",100,-180,180);
 462   fTrackDirLambda2All=new TH1F("fTrackDirLambda2All","fTrackDirLambda2All",100,0,180);
 463   fTrackDirAlphaAll=new TH1F("fTrackDirAlphaAll","fTrackDirAlphaAll",100,-180,180);
 464
 465   fTrackDirAll=new TH2F("fTrackDirAll","Hist with trakcs directions",100,-180,180,100,-180,180);
 466   fTrackDir2All=new TH2F("fTrackDir2All","Hist with trakcs directions",100,-180,180,100,-180,180);
 467  fTrackDirXZAll=new TH2F("fTrackDirXZAll","Hist with trakcs directions from XZ ",100,-3.,3.,100,-3.,3.);
 468
 469   fTrackDir=new TH2F*[fnHist];
 470
 471   Float_t **binningZPhi=CheckSingleLayer(volIDs);
 472   fvolidsToBin=new Int_t*[fnPhi*fnZ];
 473
 474   Float_t *binningphi=binningZPhi[0];
 475   Float_t *binningz=binningZPhi[1];
 476   Bool_t binning=SetBinning(volIDs,binningphi,binningz);
 477
 478   if(binning){
 479     fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
 480     fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",fnPhi,binningphi,fnZ,binningz);
 481     fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
 482     fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",fnPhi,binningphi,fnZ,binningz);
 483     fhEmpty=new TH2F("fhEmpty","Hist for getting binning",fnPhi,binningphi,fnZ,binningz);
 484     fVolNTracks->SetXTitle("Volume #phi");
 485     fVolNTracks->SetYTitle("Volume z ");
 486     fhistVolNptsUsed->SetXTitle("Volume #phi");
 487     fhistVolNptsUsed->SetYTitle("Volume z ");
 488     fhistVolUsed->SetXTitle("Volume #phi");
 489     fhistVolUsed->SetYTitle("Volume z ");
 490     fSigmaVolZ->SetXTitle("Volume #phi");
 491     fSigmaVolZ->SetYTitle("Volume z ");
 492   }
 493   else{
 494     fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
 495     fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",50,-3.2,3.2,100,-80,80);
 496     fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
 497     fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",50,-3.2,3.2,100,-80,80);
 498     fhEmpty=new TH2F("fhEmpty","Hist for getting binning",50,-3.2,3.2,100,-80,80);
 499     fVolNTracks->SetXTitle("Volume #phi");
 500     fVolNTracks->SetYTitle("Volume z ");
 501     fhistVolNptsUsed->SetXTitle("Volume #phi");
 502     fhistVolNptsUsed->SetYTitle("Volume z ");
 503     fhistVolUsed->SetXTitle("Volume #phi");
 504     fhistVolUsed->SetYTitle("Volume z ");
 505     fSigmaVolZ->SetXTitle("Volume #phi");
 506     fSigmaVolZ->SetYTitle("Volume z ");
 507   }
 508   fpTrackVolIDs=new TArrayI(fnHist);
 509   fVolUsed=new TArrayI*[fnHist];
 510   fVolVolids=new TArrayI*[fnHist];
 511   fLastVolVolid=new Int_t[fnHist];
 512
 513   for (Int_t nhist=0;nhist<fnHist;nhist++){
 514     fpTrackVolIDs->AddAt(volIDs->At(nhist),nhist);
 515     aux=histnameRPHI;
 516     aux+=volIDs->At(nhist);
 517     fVolResHistRPHI[nhist]=new TH1F("histname","histname",200,-0.02,0.02);
 518     fVolResHistRPHI[nhist]->SetName(aux.Data());
 519     fVolResHistRPHI[nhist]->SetTitle(aux.Data());
 520
 521     aux=histnameZ;
 522     aux+=volIDs->At(nhist);
 523     fResHistZ[nhist]=new TH1F("histname","histname",400,-0.08,0.08);
 524     fResHistZ[nhist]->SetName(aux.Data());
 525     fResHistZ[nhist]->SetTitle(aux.Data());
 526
 527     aux=histnamePullRPHI;
 528     aux+=volIDs->At(nhist);
 529     fPullHistRPHI[nhist]=new TH1F("histname","histname",100,-7.,7.);
 530     fPullHistRPHI[nhist]->SetName(aux.Data());
 531     fPullHistRPHI[nhist]->SetTitle(aux.Data());
 532
 533     aux=histnamePullZ;
 534     aux+=volIDs->At(nhist);
 535     fPullHistZ[nhist]=new TH1F("histname","histname",100,-7.,7.);
 536     fPullHistZ[nhist]->SetName(aux.Data());
 537     fPullHistZ[nhist]->SetTitle(aux.Data());
 538
 539     aux=histnameDirPhi;
 540     aux+=volIDs->At(nhist);
 541     fTrackDirPhi[nhist]=new TH1F("histname","histname",100,-180,180);
 542     fTrackDirPhi[nhist]->SetName(aux.Data());
 543     fTrackDirPhi[nhist]->SetTitle(aux.Data());
 544
 545     aux=histnameDirLambda;
 546     aux+=volIDs->At(nhist);
 547     fTrackDirLambda[nhist]=new TH1F("histname","histname",100,0,180);
 548     fTrackDirLambda[nhist]->SetName(aux.Data());
 549     fTrackDirLambda[nhist]->SetTitle(aux.Data());
 550
 551     aux=histnameDirLambda2;
 552     aux+=volIDs->At(nhist);
 553     fTrackDirLambda2[nhist]=new TH1F("histname","histname",100,0,180);
 554     fTrackDirLambda2[nhist]->SetName(aux.Data());
 555     fTrackDirLambda2[nhist]->SetTitle(aux.Data());
 556
 557     aux=histnameDirAlpha;
 558     aux+=volIDs->At(nhist);
 559     fTrackDirAlpha[nhist]=new TH1F("histname","histname",100,-180,180);
 560     fTrackDirAlpha[nhist]->SetName(aux.Data());
 561     fTrackDirAlpha[nhist]->SetTitle(aux.Data());
 562
 563     aux=histnameDir;
 564     aux+=volIDs->At(nhist);
 565     fTrackDir[nhist]=new TH2F("histname","histname",100,-90.,90.,100,-180.,180.);
 566     fTrackDir[nhist]->SetName(aux.Data());
 567     fTrackDir[nhist]->SetTitle(aux.Data());
 568
 569     aux=histnameGlob;
 570     aux+=volIDs->At(nhist);
 571     fResHistGlob[nhist]=new TH1F("histname","histname",400,-0.08,0.08);
 572     fResHistGlob[nhist]->SetName(aux.Data());
 573     fResHistGlob[nhist]->SetTitle(aux.Data());
 574
 575     aux=histnameCorrVol;
 576     aux+=volIDs->At(nhist);
 577     fhistCorrVol[nhist]=new TH2F("histname","histname",50,-3.2,3.2,100,-80,80);
 578
 579
 580     fhistCorrVol[nhist]->SetName(aux.Data());
 581     fhistCorrVol[nhist]->SetTitle(aux.Data());
 582     fhistCorrVol[nhist]->SetXTitle("Volume #varphi");
 583     fhistCorrVol[nhist]->SetYTitle("Volume z ");
 584     fVolVolids[nhist]=new TArrayI(100);
 585     fVolUsed[nhist]=new TArrayI(1000);
 586     fLastVolVolid[nhist]=0;
 587
 588   }
 589   fWriteHist=kFALSE;
 590
 591   return;
 592 }
 593
 594 //____________________________________________________________________________
 595 void AliITSResidualsAnalysis::ListVolUsed(TTree *pointsTree,TArrayI ***arrayIndex,Int_t **lastIndex)
 596 {
 597   //
 598   // This is copied from AliAlignmentClass::LoadPoints() method
 599   //
 600
 601   Int_t volIDalignable,volIDpoint,iModule;
 602   AliTrackPoint p;
 603   AliTrackPointArray* array = 0;
 604   pointsTree->SetBranchAddress("SP", &array);
 605
 606
 607   for(Int_t ivol=0;ivol<fnHist;ivol++){
 608     Int_t lastused=0;
 609     volIDalignable=fpTrackVolIDs->At(ivol);
 610     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volIDalignable,iModule);
 611
 612     Int_t nArraysId = lastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 613     printf("volume %d (Layer %d, Modulo %d) , numero di elementi per volume %d \n",volIDalignable,iLayer,iModule,nArraysId);
 614     TArrayI *index = arrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 615     for (Int_t iArrayId = 0;iArrayId < nArraysId; iArrayId++) {
 616
 617       // Get tree entry
 618       Int_t entry = (*index)[iArrayId];
 619
 620       pointsTree->GetEvent(entry);
 621       if (!array) {
 622         AliWarning("Wrong space point array index!");
 623         continue;
 624       }
 625
 626       // Get the space-point array
 627       Int_t modnum,nPoints = array->GetNPoints();
 628
 629       for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
 630         array->GetPoint(p,iPoint);
 631
 632         AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
 633         // check if the layer id is valid
 634         if ((layer < AliGeomManager::kFirstLayer) ||
 635             (layer >= AliGeomManager::kLastLayer)) {
 636           AliError(Form("Layer index is invalid: %d (%d -> %d) !",
 637                         layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
 638           continue;
 639         }
 640         if ((modnum >= AliGeomManager::LayerSize(layer)) ||
 641             (modnum < 0)) {
 642           AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
 643                         layer,modnum,AliGeomManager::LayerSize(layer)));
 644           continue;
 645         }
 646         if (layer > AliGeomManager::kSSD2) continue; // ITS only
 647
 648         volIDpoint=(Int_t)p.GetVolumeID();
 649         if (volIDpoint==volIDalignable)continue;
 650         Int_t size = fVolVolids[ivol]->GetSize();
 651         // If needed allocate new size
 652         if (fLastVolVolid[ivol]>=size){// Warning: fLAST[NHIST] is useless
 653           fVolVolids[ivol]->Set(size + 1000);
 654         }
 655         fVolVolids[ivol]->AddAt(volIDpoint,fLastVolVolid[ivol]);
 656         fLastVolVolid[ivol]++;
 657         Bool_t usedVol=kFALSE;
 658         for(Int_t used=0;used<lastused;used++){
 659           if(fVolUsed[ivol]->At(used)==volIDpoint){
 660             usedVol=kTRUE;
 661             break;
 662           }
 663         }
 664         if (!usedVol){
 665           size = fVolUsed[ivol]->GetSize();
 666           // If needed allocate new size
 667           if (lastused>= size){
 668             fVolUsed[ivol]->Set(size + 1000);
 669           }
 670           fVolUsed[ivol]->AddAt(volIDpoint,lastused);
 671           lastused++;
 672         }
 673
 674         FillVolumeCorrelationHists(ivol,volIDalignable,volIDpoint,usedVol);
 675       }
 676     }
 677   }
 678   fWriteHist=kTRUE;
 679   return;
 680 }
 681
 682 //____________________________________________________________________________
 683 void AliITSResidualsAnalysis::FillVolumeCorrelationHists(Int_t ivol,Int_t volIDalignable,Int_t volIDpoint,Bool_t usedVol) const
 684 {
 685   //
 686   // Fill the histograms with the correlations between volumes
 687   //
 688
 689   if(!gGeoManager)AliGeomManager::LoadGeometry(GetFileNameGeometry());
 690   Double_t *transGlobal,radius,phi;
 691   const char *symname,*volpath;
 692   TGeoPNEntry *pne;
 693   TGeoPhysicalNode *pn;
 694   TGeoHMatrix *globMatrix;
 695
 696   symname = AliGeomManager::SymName(volIDalignable);
 697   pne = gGeoManager->GetAlignableEntry(symname);
 698   volpath=pne->GetTitle();
 699   pn=gGeoManager->MakePhysicalNode(volpath);
 700   globMatrix=pn->GetMatrix();
 701
 702   transGlobal=globMatrix->GetTranslation();
 703   radius=TMath::Sqrt(transGlobal[0]*transGlobal[0]+transGlobal[1]*transGlobal[1]);
 704   phi=TMath::ATan2(transGlobal[1],transGlobal[0]);
 705   fhistVolNptsUsed->Fill(phi,transGlobal[2]);
 706   if(!usedVol)fhistVolUsed->Fill(phi,transGlobal[2]);
 707
 708   symname = AliGeomManager::SymName(volIDpoint);
 709   pne = gGeoManager->GetAlignableEntry(symname);
 710   volpath=pne->GetTitle();
 711   pn=gGeoManager->MakePhysicalNode(volpath);
 712   globMatrix=pn->GetMatrix();
 713
 714   transGlobal=globMatrix->GetTranslation();
 715   radius=TMath::Sqrt(transGlobal[0]*transGlobal[0]+transGlobal[1]*transGlobal[1]);
 716   phi=TMath::ATan2(transGlobal[1],transGlobal[0]);
 717
 718   fhistCorrVol[ivol]->Fill(phi,transGlobal[2]);
 719
 720   return;
 721 }
 722
 723
 724 //____________________________________________________________________________
 725 void AliITSResidualsAnalysis::FillResHists(AliTrackPointArray *points,AliTrackPointArray *pTrack) const
 726 {
 727   //
 728   // Method that fills the histograms with the residuals
 729   //
 730
 731   Int_t volIDpoint;
 732   Float_t xyz[3],xyz2[3];
 733   const Float_t *cov,*cov2;
 734   Float_t resRPHI,resGlob,resZ;
 735   Double_t pullz,pullrphi,sign;
 736   Double_t phi,lambda,lambda2,alpha,xovery,zovery;
 737   AliTrackPoint p,pTr;
 738   for(Int_t ipoint=0;ipoint<points->GetNPoints();ipoint++){
 739     points->GetPoint(p,ipoint);
 740     volIDpoint=(Int_t)p.GetVolumeID();
 741     p.GetXYZ(xyz);
 742     cov=p.GetCov();
 743     pTrack->GetPoint(pTr,ipoint);
 744     GetTrackDirClusterCov(&pTr,phi,lambda,lambda2,alpha,xovery,zovery);
 745     pTr.GetXYZ(xyz2);
 746     cov2=pTr.GetCov();
 747     resRPHI=TMath::Sqrt((xyz2[0]-xyz[0])*(xyz2[0]-xyz[0])+(xyz2[1]-xyz[1])*(xyz2[1]-xyz[1]));
 748     //resRPHI is always positive value
 749     sign=TMath::ATan2(xyz2[1],xyz2[0])-TMath::ATan2(xyz[1],xyz[0]);
 750     if(sign!=0.){
 751       sign=sign/TMath::Abs(sign);
 752       resRPHI=resRPHI*sign;
 753       pullrphi=sign*resRPHI*resRPHI/TMath::Sqrt((xyz2[0]-xyz[0])*(xyz2[0]-xyz[0])*(cov2[0]/100000000.+cov[0])+(xyz2[1]-xyz[1])*(xyz2[1]-xyz[1])*(cov2[3]/100000000.+cov[3]));
 754     }
 755     else{
 756       pullrphi=0.;
 757       resRPHI=0.;
 758     }
 759
 760     resZ=xyz2[2]-xyz[2];
 761     pullz=resZ/(TMath::Sqrt(cov2[5])/10000.);
 762     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]));
 763     for(Int_t ivolIDs=0;ivolIDs<fpTrackVolIDs->GetSize();ivolIDs++){
 764       if(volIDpoint==fpTrackVolIDs->At(ivolIDs)){
 765         fVolResHistRPHI[ivolIDs]->Fill(resRPHI);
 766         fResHistZ[ivolIDs]->Fill(resZ);
 767         fResHistGlob[ivolIDs]->Fill(resGlob);
 768
 769
 770         fTrackDirPhi[ivolIDs]->Fill(phi);
 771         fTrackDirLambda[ivolIDs]->Fill(lambda);
 772         fTrackDirLambda2[ivolIDs]->Fill(lambda2);
 773         fTrackDirAlpha[ivolIDs]->Fill(alpha);
 774
 775         fTrackDirPhiAll->Fill(phi);
 776         fTrackDirLambdaAll->Fill(lambda);
 777         fTrackDirLambda2All->Fill(lambda2);
 778         fTrackDirAlphaAll->Fill(alpha);
 779
 780         fTrackDirAll->Fill(lambda,alpha);
 781         fTrackDir2All->Fill(lambda2,phi);
 782         fTrackDirXZAll->Fill(xovery,zovery);
 783         fTrackDir[ivolIDs]->Fill(lambda,alpha);
 784
 785         fPullHistRPHI[ivolIDs]->Fill(pullrphi);
 786         fPullHistZ[ivolIDs]->Fill(pullz);
 787
 788         if(fsingleLayer){
 789           Int_t binz,binphi;
 790           Float_t globalPhi,globalZ;
 791           if(kTRUE||(fvolidsToBin[ivolIDs][0]!=volIDpoint)){
 792             binphi=GetBinPhiZ((Int_t)volIDpoint,&binz);
 793           }
 794           else{//this in the case of alignment of one entire layer (fnHIst=layersize) may reduce iterations: remind of that fsingleLayer->fnHista<layerSize
 795             binphi=fvolidsToBin[ivolIDs][1];
 796             binz=fvolidsToBin[ivolIDs][2];
 797           }
 798           globalPhi=fCoordToBinTable[binphi][binz][0];
 799           globalZ=fCoordToBinTable[binphi][binz][1];
 800
 801           fVolNTracks->Fill(globalPhi,globalZ);
 802         }
 803         else fVolNTracks->Fill(TMath::ATan2(xyz[1],xyz[0]),xyz[2]);
 804       }
 805     }
 806   }
 807 }
 808
 809
 810 //____________________________________________________________________________
 811 Bool_t AliITSResidualsAnalysis::AnalyzeHists(Int_t minNpoints) const
 812 {
 813   //
 814   // Saves the histograms into a tree and saves the tree into a file
 815   //
 816
 817   TString outname = "ResidualsAnalysisTree.root";
 818   TFile *hFile=new TFile(outname.Data(),"RECREATE","The Files containing the TREE with Align. Vol. hists nd Prop.");
 819   TTree *analysisTree=new TTree("analysisTree","Tree whith residuals analysis data for alignable volumes");
 820   static histProperties_t histRPHIprop,histZprop,histGlobprop;
 821   Int_t volID;
 822
 823   TF1 *gauss;
 824   TH1F *histRPHI,*histZ,*histGlob,*histPullRPHI,*histPullZ,*hTrackDirPhi,*hTrackDirLambda,*hTrackDirLambda2,*hTrackDirAlpha;
 825
 826   TH2F *histCorrVol,*hTrackDir;
 827
 828   histRPHI=new TH1F();
 829   histZ=new TH1F();
 830   histPullRPHI=new TH1F();
 831   histPullZ=new TH1F();
 832   hTrackDirPhi=new TH1F();
 833   hTrackDirLambda=new TH1F();
 834   hTrackDirLambda2=new TH1F();
 835   hTrackDirAlpha=new TH1F();
 836   hTrackDir=new TH2F();
 837   histGlob=new TH1F();
 838   histCorrVol=new TH2F();
 839   Float_t globalPhi,globalZ;
 840   Double_t rms;
 841   Int_t nHistAnalyzed=0,entries;
 842   analysisTree->Branch("volID",&volID,"volID/I");
 843   if(fsingleLayer){
 844     analysisTree->Branch("globalPhi",&globalPhi,"globalPhi/F");
 845     analysisTree->Branch("globalZ",&globalZ,"globalZ/F");
 846   }
 847   analysisTree->Branch("histRPHI","TH1F",&histRPHI,128000,0);
 848   analysisTree->Branch("histPullRPHI","TH1F",&histPullRPHI,128000,0);
 849
 850   analysisTree->Branch("histRPHIprop",&histRPHIprop,"nentries/I:rms/F:meanFit:errmeanFit:sigmaFit");
 851   analysisTree->Branch("histZ","TH1F",&histZ,128000,0);
 852   analysisTree->Branch("histPullZ","TH1F",&histPullZ,128000,0);
 853   analysisTree->Branch("hTrackDirPhi","TH1F",&hTrackDirPhi,128000,0);
 854   analysisTree->Branch("hTrackDirLambda","TH1F",&hTrackDirLambda,128000,0);
 855   analysisTree->Branch("hTrackDirLambda2","TH1F",&hTrackDirLambda2,128000,0);
 856   analysisTree->Branch("hTrackDirAlpha","TH1F",&hTrackDirAlpha,128000,0);
 857   analysisTree->Branch("hTrackDir","TH2F",&hTrackDir,128000,0);
 858
 859   analysisTree->Branch("histZprop",&histZprop,"nentries/I:rms/F:meanFit:errmeanFit:sigmaFit");
 860   analysisTree->Branch("histGlob","TH1F",&histGlob,128000,0);
 861   analysisTree->Branch("histGlobprop",&histGlobprop,"nentries/I:rms/F:meanFit:errmeanFit:sigmaFit");
 862   if(fWriteHist){
 863     analysisTree->Branch("histCorrVol","TH2F",&histCorrVol,128000,0);
 864   }
 865
 866   for(Int_t j=0;j<fnHist;j++){
 867     volID=fpTrackVolIDs->At(j);
 868     if((entries=(fResHistGlob[j]->GetEntries())>=minNpoints)||fsingleLayer){
 869       nHistAnalyzed++;
 870       //histRPHI
 871       histRPHI=fVolResHistRPHI[j];
 872       histPullRPHI=fPullHistRPHI[j];
 873       histRPHIprop.nentries=(Int_t)fVolResHistRPHI[j]->GetEntries();
 874       rms=fVolResHistRPHI[j]->GetRMS();
 875       histRPHIprop.rms=rms;
 876       gauss=new TF1("gauss","gaus",-3*rms,3*rms);
 877       fVolResHistRPHI[j]->Fit("gauss","RN");
 878       histRPHIprop.meanFit=gauss->GetParameter(1);
 879       histRPHIprop.errmeanFit=gauss->GetParError(1);
 880       histRPHIprop.sigmaFit=gauss->GetParameter(2);
 881       //histZ
 882       histZ=fResHistZ[j];
 883       histPullZ=fPullHistZ[j];
 884       histZprop.nentries=(Int_t)fResHistZ[j]->GetEntries();
 885       rms=fResHistZ[j]->GetRMS();
 886       histZprop.rms=rms;
 887       gauss=new TF1("gauss","gaus",-3*rms,3*rms);
 888       fResHistZ[j]->Fit("gauss","RN");
 889       histZprop.meanFit=gauss->GetParameter(1);
 890       histZprop.errmeanFit=gauss->GetParError(1);
 891       histZprop.sigmaFit=gauss->GetParameter(2);
 892       //histGlob
 893       histGlob=fResHistGlob[j];
 894       histGlobprop.nentries=(Int_t)fResHistGlob[j]->GetEntries();
 895       rms=fResHistGlob[j]->GetRMS();
 896       histGlobprop.rms=rms;
 897       gauss=new TF1("gauss","gaus",-3*rms,3*rms);
 898       fResHistGlob[j]->Fit("gauss","RN");
 899       histGlobprop.meanFit=gauss->GetParameter(1);
 900       histGlobprop.errmeanFit=gauss->GetParError(1);
 901       histGlobprop.sigmaFit=gauss->GetParameter(2);
 902
 903       //histTrackDir
 904       hTrackDirPhi=fTrackDirPhi[j];
 905       hTrackDirLambda=fTrackDirLambda[j];
 906       hTrackDirLambda2=fTrackDirLambda2[j];
 907       hTrackDirAlpha=fTrackDirAlpha[j];
 908       hTrackDir=fTrackDir[j];
 909
 910       if(fsingleLayer){
 911         Int_t binz,binphi;
 912         if (fvolidsToBin[j][0]!=volID)binphi=GetBinPhiZ((Int_t)volID,&binz);
 913         else{
 914           binphi=fvolidsToBin[j][1];
 915           binz=fvolidsToBin[j][2];
 916         }
 917         globalPhi=fCoordToBinTable[binphi][binz][0];
 918         globalZ=fCoordToBinTable[binphi][binz][1];
 919
 920
 921         histCorrVol=fhistCorrVol[j];
 922         fSigmaVolZ->SetBinContent(binphi+1,binz+1,histRPHIprop.sigmaFit);//+1 is for underflow
 923       }
 924       analysisTree->Fill();
 925     }
 926     else continue;
 927
 928   }
 929   if(nHistAnalyzed>0){
 930     analysisTree->Print();
 931     fVolNTracks->Write();
 932     hFile->Write();
 933     fhEmpty->Write();
 934     if(fWriteHist){
 935       fhistVolUsed->Draw();
 936       fSigmaVolZ->Draw();
 937       fSigmaVolZ->Write();
 938       fhistVolUsed->Write();
 939       fTrackDirPhiAll->Write();
 940       fTrackDirLambdaAll->Write();
 941       fTrackDirLambda2All->Write();
 942       fTrackDirAlphaAll->Write();
 943       fTrackDirAll->Write();
 944       fTrackDir2All->Write();
 945       fTrackDirXZAll->Write();
 946       fhistVolNptsUsed->Write();
 947       hFile->Close();
 948     }
 949     return kTRUE;
 950   }
 951   else {
 952     delete analysisTree;
 953     delete hFile;
 954     return kFALSE;}
 955 }
 956
 957
 958 //____________________________________________________________________________
 959 void AliITSResidualsAnalysis::DrawHists() const
 960 {
 961   //
 962   // Draws the histograms of the residuals and of the number of tracks
 963   //
 964
 965   TString cname;
 966   for(Int_t canv=0;canv<fnHist;canv++){
 967     cname="canv Residuals";
 968     cname+=canv;
 969     TCanvas *c=new TCanvas(cname.Data(),cname.Data(),700,700);
 970     c->Divide(3,1);
 971     c->cd(1);
 972     fVolResHistRPHI[canv]->Draw();
 973     c->cd(2);
 974     fResHistZ[canv]->Draw();
 975     c->cd(3);
 976     fResHistGlob[canv]->Draw();
 977   }
 978   cname="canv NVolTracks";
 979
 980   TCanvas *c2=new TCanvas(cname.Data(),cname.Data(),700,700);
 981   c2->cd();
 982   fVolNTracks->Draw();
 983
 984   return;
 985 }
 986
 987 //____________________________________________________________________________
 988 Float_t** AliITSResidualsAnalysis::CheckSingleLayer(const TArrayI *volids)
 989 {
 990   //
 991   // Checks if volumes array is a single (ITS) layer or not
 992   //
 993
 994   Float_t **binningzphi=new Float_t*[2];
 995   Int_t iModule;
 996   AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volids->At(0),iModule);
 997   //Check that one single Layer is going to be aligned
 998   for(Int_t nvol=0;nvol<volids->GetSize();nvol++){
 999     if(iLayer != AliGeomManager::VolUIDToLayer((UShort_t)volids->At(nvol),iModule)){
1000       printf("Wrong Layer! \n %d , %d , %d ,%d \n",(Int_t)AliGeomManager::VolUIDToLayer((UShort_t)volids->At(nvol),iModule),nvol,volids->GetSize(),iModule);
1001       fsingleLayer=kFALSE;
1002       return 0;
1003     }
1004   }
1005
1006   //Bool_t used=kFALSE;
1007   switch (iLayer) {
1008   case AliGeomManager::kSPD1:{
1009     fnPhi=fkPhiSPD1;//fkPhiSPD1;
1010     fnZ=fkZSPD1;//nZSPD1;
1011     binningzphi[0]=new Float_t[fkPhiSPD1+1];
1012     binningzphi[1]=new Float_t[fkZSPD1+1];
1013     fCoordToBinTable=new Double_t**[fkPhiSPD1];
1014     for(Int_t j=0;j<fnPhi;j++){
1015       fCoordToBinTable[j]=new Double_t*[fkZSPD1];
1016     }
1017   }; break;
1018   case AliGeomManager::kSPD2:{
1019     fnPhi=fkPhiSPD2;//fkPhiSPD1;
1020     fnZ=fkZSPD2;//nZSPD1;
1021     binningzphi[0]=new Float_t[fkPhiSPD2+1];
1022     binningzphi[1]=new Float_t[fkZSPD2+1];
1023     fCoordToBinTable=new Double_t**[fkPhiSPD2];
1024     for(Int_t j=0;j<fnPhi;j++){
1025       fCoordToBinTable[j]=new Double_t*[fkZSPD2];
1026     }
1027
1028   }; break; case AliGeomManager::kSDD1:{
1029     fnPhi=fkPhiSDD1;//fkPhiSPD1;
1030     fnZ=fkZSDD1;//nZSPD1;
1031     binningzphi[0]=new Float_t[fkPhiSDD1+1];
1032     binningzphi[1]=new Float_t[fkZSDD1+1];
1033     fCoordToBinTable=new Double_t**[fkPhiSDD1];
1034     for(Int_t j=0;j<fnPhi;j++){
1035       fCoordToBinTable[j]=new Double_t*[fkZSDD1];
1036     }
1037   }; break; case AliGeomManager::kSDD2:{
1038     fnPhi=fkPhiSDD2;//fkPhiSPD1;
1039     fnZ=fkZSDD2;//nZSPD1;
1040     binningzphi[0]=new Float_t[fkPhiSDD2+1];
1041     binningzphi[1]=new Float_t[fkZSDD2+1];
1042     fCoordToBinTable=new Double_t**[fkPhiSDD2];
1043     for(Int_t j=0;j<fnPhi;j++){
1044       fCoordToBinTable[j]=new Double_t*[fkZSDD2];
1045     }
1046   }; break; case AliGeomManager::kSSD1:{
1047     fnPhi=fkPhiSSD1;//fkPhiSPD1;
1048     fnZ=fkZSSD1;//nZSPD1;
1049     binningzphi[0]=new Float_t[fkPhiSSD1+1];
1050     binningzphi[1]=new Float_t[fkZSSD1+1];
1051     fCoordToBinTable=new Double_t**[fkPhiSSD1];
1052     for(Int_t j=0;j<fnPhi;j++){
1053       fCoordToBinTable[j]=new Double_t*[fkZSSD1];
1054     }
1055   }; break; case AliGeomManager::kSSD2:{
1056     fnPhi=fkPhiSSD2;//fkPhiSPD1;
1057     fnZ=fkZSSD2;//nZSPD1;
1058     binningzphi[0]=new Float_t[fkPhiSSD2+1];
1059     binningzphi[1]=new Float_t[fkZSSD2+1];
1060     fCoordToBinTable=new Double_t**[fkPhiSSD2];
1061     for(Int_t j=0;j<fnPhi;j++){
1062       fCoordToBinTable[j]=new Double_t*[fkZSSD2];
1063     }
1064   }; break;
1065
1066   default:{
1067     printf("Wrong Layer Label! \n");
1068     fsingleLayer=kFALSE;
1069     return 0x0;
1070   }
1071   }
1072   fsingleLayer=kTRUE;
1073   return binningzphi;
1074 }
1075
1076 //____________________________________________________________________________
1077 Bool_t AliITSResidualsAnalysis::SetBinning(const TArrayI *volids,Float_t *phiBin,Float_t *zBin)
1078 {
1079   //
1080   // Sets the correct binning
1081   //
1082
1083   if(!fsingleLayer)return kFALSE;
1084   const char *volpath,*symname;
1085   Int_t iModule;
1086   Int_t *orderArrayPhi,*orderArrayZ;
1087   UShort_t volID;
1088   Double_t *phiArray,*zArray,*transGlobal,*phiArrayOrdered,*zArrayOrdered;
1089   Double_t lastPhimin=-10;
1090   Double_t lastZmin=-99;
1091   Int_t ***orderPhiZ;
1092   TGeoPNEntry *pne;
1093   TGeoPhysicalNode *pn;
1094   TGeoHMatrix *globMatrix;
1095
1096   Bool_t used=kFALSE;
1097
1098   orderPhiZ=new Int_t**[fnPhi];
1099   phiArray=new Double_t[fnPhi];//phiBin[nModulesPhi+1];
1100   zArray=new Double_t[fnZ];//zBin[nModulesZ+1];
1101   phiArrayOrdered=new Double_t[fnPhi];
1102   zArrayOrdered=new Double_t[fnZ];
1103   orderArrayPhi=new Int_t[fnPhi];
1104   orderArrayZ=new Int_t[fnZ];
1105   for(Int_t k=0;k<fnZ;k++){
1106     orderArrayZ[k]=0;
1107     zArray[k]=-1000;
1108   }
1109   for(Int_t k=0;k<fnPhi;k++){
1110     orderArrayPhi[k]=0;
1111     phiArray[k]=-10;
1112     orderPhiZ[k]=new Int_t*[fnZ];
1113   }
1114
1115
1116   AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volids->At(0),iModule);
1117
1118   Int_t lastPhi=0,lastZ=0;
1119   for(iModule=0;iModule<AliGeomManager::LayerSize(iLayer);iModule++){
1120     fvolidsToBin[iModule]=new Int_t[3];
1121     volID=AliGeomManager::LayerToVolUID(iLayer,iModule);
1122     fvolidsToBin[iModule][0]=volID;
1123     symname = AliGeomManager::SymName(volID);
1124     pne = gGeoManager->GetAlignableEntry(symname);
1125     volpath=pne->GetTitle();
1126     pn=gGeoManager->MakePhysicalNode(volpath);
1127     globMatrix=pn->GetMatrix();
1128     transGlobal=globMatrix->GetTranslation();
1129
1130     for(Int_t j=0;j<lastPhi;j++){
1131       used=kFALSE;
1132       if(TMath::Abs(phiArray[j]-TMath::ATan2(transGlobal[1],transGlobal[0]))<2*TMath::Pi()/(10*fnPhi)){//10 is a safety factor but....
1133         fvolidsToBin[iModule][1]=j;
1134         used=kTRUE;
1135         break;
1136       }
1137     }
1138     if(!used){
1139       phiArray[lastPhi]=TMath::ATan2(transGlobal[1],transGlobal[0]);
1140       fvolidsToBin[iModule][1]=lastPhi;
1141       if(phiArray[lastPhi]<lastPhimin)lastPhimin=phiArray[lastPhi];
1142       lastPhi++;
1143       if(lastPhi>fnPhi){
1144         printf("Wrong Phi! \n");
1145         return kFALSE;}
1146     }
1147
1148     for(Int_t j=0;j<lastZ;j++){
1149       used=kFALSE;
1150       if(TMath::Abs(zArray[j]-transGlobal[2])<0.1){
1151         fvolidsToBin[iModule][2]=j;
1152         used=kTRUE;
1153         break;
1154       }
1155     }
1156     if(!used){
1157       fvolidsToBin[iModule][2]=lastZ;
1158       zArray[lastZ]=transGlobal[2];
1159       if(zArray[lastZ]<lastZmin)lastZmin=zArray[lastZ];
1160       lastZ++;
1161       if(lastZ>fnZ){
1162         printf("Wrong Z! \n");
1163         return kFALSE;}
1164     }
1165   }
1166
1167
1168   //ORDERING THE ARRAY OF PHI AND Z VALUES
1169   for(Int_t order=0;order<fnPhi;order++){
1170     for(Int_t j=0;j<fnPhi;j++){
1171       if((j!=order)&&(phiArray[order]>phiArray[j]))orderArrayPhi[order]++;
1172     }
1173   }
1174
1175   for(Int_t order=0;order<fnPhi;order++){
1176     for(Int_t j=0;j<fnPhi;j++){
1177       if(orderArrayPhi[j]==order){
1178         phiArrayOrdered[order]=phiArray[j];
1179         break;
1180         }
1181     }
1182   }
1183
1184
1185   for(Int_t order=0;order<fnZ;order++){
1186     for(Int_t j=0;j<fnZ;j++){
1187       if((j!=order)&&(zArray[order]>zArray[j]))orderArrayZ[order]++;
1188     }
1189   }
1190
1191
1192   for(Int_t order=0;order<fnZ;order++){
1193     for(Int_t j=0;j<fnZ;j++){
1194       if(orderArrayZ[j]==order){
1195         zArrayOrdered[order]=zArray[j];
1196         break;
1197       }
1198     }
1199   }
1200
1201
1202   //Filling the  fCoordToBinTable
1203   for(Int_t j=0;j<fnPhi;j++){
1204     for(Int_t i=0;i<fnZ;i++){
1205       orderPhiZ[j][i]=new Int_t[2];
1206       orderPhiZ[j][i][0]=orderArrayPhi[j];
1207       orderPhiZ[j][i][1]=orderArrayZ[i];
1208     }
1209   }
1210
1211
1212   for(Int_t j=0;j<fnPhi;j++){
1213     for(Int_t i=0;i<fnZ;i++){
1214       fCoordToBinTable[j][i]=new Double_t[2];
1215       fCoordToBinTable[j][i][0]=phiArrayOrdered[j];
1216       fCoordToBinTable[j][i][1]=zArrayOrdered[i];
1217       printf("ecco (binphi,binz)= %d, %d e (phi,z)=%f,%f \n",j,i,fCoordToBinTable[j][i][0],fCoordToBinTable[j][i][1]);
1218     }
1219   }
1220   Int_t istar,jstar;
1221   for(iModule=0;iModule<fnPhi*fnZ;iModule++){
1222     istar=fvolidsToBin[iModule][1];
1223     jstar=fvolidsToBin[iModule][2];
1224     fvolidsToBin[iModule][1]=orderPhiZ[istar][jstar][0];
1225     fvolidsToBin[iModule][2]=orderPhiZ[istar][jstar][1];
1226   }
1227
1228
1229   //now constructing the binning
1230   for(Int_t iModPhi=0;iModPhi<fnPhi-1;iModPhi++){
1231     phiBin[iModPhi+1]=(Float_t)phiArrayOrdered[iModPhi]+0.5*(phiArrayOrdered[iModPhi+1]-phiArrayOrdered[iModPhi]);
1232   }
1233
1234   phiBin[0]=(Float_t)phiArrayOrdered[0]-(phiArrayOrdered[1]-phiArrayOrdered[0])/2.;
1235
1236   phiBin[fnPhi]=(Float_t)phiArrayOrdered[fnPhi-1]+(phiArrayOrdered[fnPhi-1]-phiArrayOrdered[fnPhi-2])/2.;
1237   for(Int_t iModPhi=0;iModPhi<fnPhi+1;iModPhi++){
1238     printf("Mean Phi mod %d su %d:  %f \n",iModPhi+1,fnPhi,phiBin[iModPhi]);
1239   }
1240
1241   for(Int_t iModZ=0;iModZ<fnZ-1;iModZ++){
1242     zBin[iModZ+1]=(Float_t)zArrayOrdered[iModZ]+0.5*(zArrayOrdered[iModZ+1]-zArrayOrdered[iModZ]);
1243   }
1244   zBin[0]=(Float_t)zArrayOrdered[0]-0.5*(zArrayOrdered[1]-zArrayOrdered[0]);
1245   zBin[fnZ]=(Float_t)zArrayOrdered[fnZ-1]+0.5*(zArrayOrdered[1]-zArrayOrdered[0]);
1246
1247
1248   for(Int_t iModPhi=0;iModPhi<fnZ+1;iModPhi++){
1249      printf("Mean Z mod %d su %d:  %f \n",iModPhi+1,fnZ,zBin[iModPhi]);
1250   }
1251   return kTRUE;
1252 }
1253
1254
1255 //____________________________________________________________________________
1256 Int_t AliITSResidualsAnalysis::GetBinPhiZ(const Int_t volID,Int_t *binz) const
1257 {
1258   //
1259   // Returns the correct Phi-Z bin
1260   //
1261
1262   if (!fsingleLayer){
1263     printf("No Single Layer reAlignment! \n");
1264     return 100;
1265   }
1266
1267   for(Int_t j=0;j<fnPhi*fnZ;j++){
1268     if(j==fnZ*fnPhi){
1269       printf("Wrong volume UID introduced! fnHist: %d  volID: %d iter: %d \n",fnHist,volID,j);
1270       return 100;
1271     }
1272     if(fvolidsToBin[j][0]==volID){
1273
1274       *binz=fvolidsToBin[j][2];
1275       return fvolidsToBin[j][1];
1276     }
1277   }
1278
1279   return 100;
1280
1281 }
1282
1283 //____________________________________________________________________________
1284 TArrayI* AliITSResidualsAnalysis::GetSingleLayerVolids(Int_t layer) const
1285 {
1286   //
1287   // Translates the layer number into a Volumes Array
1288   //
1289
1290   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1291
1292   if(layer<1 || layer>6){
1293     printf("WRONG LAYER SET! \n");
1294     return 0;
1295   }
1296   Int_t iModule,size;
1297   UShort_t volid;
1298   size = AliGeomManager::LayerSize(layer);
1299   TArrayI *volIDs = new TArrayI(size);
1300   for(iModule=0;iModule<size;iModule++){
1301     volid = AliGeomManager::LayerToVolUID(layer,iModule);
1302     volIDs->AddAt(volid,iModule);
1303   }
1304
1305   return volIDs;
1306
1307 }
1308
1309 //____________________________________________________________________________
1310 void AliITSResidualsAnalysis::GetTrackDirClusterCov(AliTrackPoint *point,Double_t &phi,Double_t &lambda,Double_t &lambda2,Double_t &alpha,Double_t &xovery,Double_t &zovery) const
1311 {
1312   //
1313   // ...
1314   //
1315
1316   TMatrixDSym cov(3);
1317   const Float_t *covvector=point->GetCov();
1318   cov(0,0)=covvector[0];
1319   cov(1,0)=cov(0,1)=covvector[1];
1320   cov(2,0)=cov(0,2)=covvector[2];
1321   cov(1,1)=covvector[3];
1322   cov(1,2)=cov(2,1)=covvector[4];
1323   cov(2,2)=covvector[5];
1324
1325   Double_t determinant=cov.Determinant();
1326   if(determinant!=0.){
1327     TMatrixD vect(3,3);
1328     TVectorD eigenvalues(3);
1329     const TMatrixDSymEigen keigen(cov);
1330     eigenvalues=keigen.GetEigenValues();
1331     vect=keigen.GetEigenVectors();
1332     Double_t mainvect[3];
1333     mainvect[0]=vect(0,0);
1334     mainvect[1]=vect(1,0);
1335     mainvect[2]=vect(2,0);
1336     if(mainvect[1]!=0.){
1337       xovery=mainvect[0]/mainvect[1];
1338       zovery=mainvect[2]/mainvect[1];
1339     }
1340     else {
1341       xovery=9999.;
1342       zovery=9999.;
1343     }
1344     if(mainvect[1]<0.){
1345       mainvect[0]=-1.*mainvect[0];
1346       mainvect[1]=-1.*mainvect[1];
1347       mainvect[2]=-1.*mainvect[2];
1348     }
1349     lambda2=TMath::ATan2(TMath::Sqrt(mainvect[0]*mainvect[0]+mainvect[2]*mainvect[2]),mainvect[1])*TMath::RadToDeg();
1350     lambda=TMath::ATan2(mainvect[2],TMath::Sqrt(mainvect[0]*mainvect[0]+mainvect[1]*mainvect[1]))*TMath::RadToDeg();
1351     phi=TMath::ATan2(mainvect[0],mainvect[2])*TMath::RadToDeg();
1352     alpha=TMath::ATan2(mainvect[1],mainvect[0])*TMath::RadToDeg();
1353   }
1354   else printf("determinant =0!, skip this point \n");
1355
1356   return;
1357 }
1358
1359 //____________________________________________________________________________
1360 void AliITSResidualsAnalysis::CalculateResiduals(const TArrayI *volids,
1361       const TArrayI *volidsfit,
1362       AliGeomManager::ELayerID layerRangeMin,
1363       AliGeomManager::ELayerID layerRangeMax,
1364       Int_t iterations,
1365       Bool_t draw)
1366 {
1367   // CalculateResiduals for a set of detector volumes.
1368   // Tracks are fitted only within
1369   // the range defined by the user
1370   // (by layerRangeMin and layerRangeMax)
1371   // or within the set of volidsfit
1372   // Repeat the procedure 'iterations' times
1373
1374
1375   Int_t nVolIds = volids->GetSize();
1376   if (nVolIds == 0) {
1377     AliError("Volume IDs array is empty!");
1378     return;
1379   }
1380
1381   // Load only the tracks with at least one
1382   // space point in the set of volume (volids)
1383
1384   //AliAlignmentTracks::SetPointsFilename(GetFileNameTrackPoints());
1385   AliAlignmentTracks::BuildIndex();
1386
1387   cout<<" Index Built!"<<endl;
1388
1389   if(draw) ListVolUsed(fPointsTree,fArrayIndex,fLastIndex);
1390
1391   AliTrackPointArray **points;
1392
1393   // Start the iterations
1394   while (iterations > 0) {
1395     Int_t nArrays = LoadPoints(volids, points);
1396     if (nArrays == 0) return;
1397
1398     AliTrackResiduals *minimizer = CreateMinimizer();
1399     minimizer->SetNTracks(nArrays);
1400     minimizer->InitAlignObj();
1401     AliTrackFitter *fitter = CreateFitter();
1402
1403     for (Int_t iArray = 0; iArray < nArrays; iArray++) {
1404       if (!points[iArray]) continue;
1405
1406
1407       fitter->SetTrackPointArray(points[iArray],kFALSE);
1408       if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
1409       AliTrackPointArray *pVolId,*pTrack;
1410
1411
1412       fitter->GetTrackResiduals(pVolId,pTrack);
1413       if(draw) FillResHists(pVolId,pTrack); // WARNING!
1414
1415       minimizer->AddTrackPointArrays(pVolId,pTrack);
1416
1417     }
1418
1419     if(minimizer->GetNFilledTracks()<=1){
1420       printf("No good tracks found: could not find parameter for volume %d (and following in volids)\n",volids->At(0));
1421       UnloadPoints(nArrays, points);
1422       return;
1423     }
1424
1425     minimizer->Minimize();
1426
1427     // Update the alignment object(s)
1428     Int_t last=0;
1429
1430     if(fRealignObjFileIsOpen)last=fClonesArray->GetLast();
1431
1432
1433     if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
1434       UShort_t volid = (*volids)[iVolId];
1435       Int_t iModule;
1436       AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
1437       AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
1438       *alignObj *= *minimizer->GetAlignObj();
1439
1440       if(iterations==1)alignObj->Print("");
1441       if(iterations==1&&fRealignObjFileIsOpen){
1442         TClonesArray &alo=*fClonesArray;
1443         new(alo[last+1+(Int_t)iVolId])AliAlignObjParams(*alignObj);
1444       }
1445
1446
1447     }
1448
1449
1450     UnloadPoints(nArrays, points);
1451
1452     iterations--;
1453
1454
1455     if(draw && iterations==0) AnalyzeHists(30);
1456
1457   }
1458
1459   return;
1460
1461 }
1462
1463
1464 //______________________________________________________________________________
1465 void AliITSResidualsAnalysis::ProcessPoints(TString minimizer,
1466       Int_t fit,
1467       AliGeomManager::ELayerID iLayerToAlign,
1468       AliGeomManager::ELayerID iLayerToExclude,
1469       TString misalignmentFile)
1470 {
1471   //
1472   // This function process the AliTrackPoints (into residuals)
1473   //
1474
1475   SetPointsFilename(GetFileNameTrackPoints());
1476   AliTrackFitter *fitter;
1477
1478   if(fit==1){
1479     fitter = new AliTrackFitterKalman();
1480   }else fitter = new AliTrackFitterRieman();
1481
1482   fitter->SetMinNPoints(4);
1483   SetTrackFitter(fitter);
1484
1485
1486   AliTrackResiduals *res;
1487
1488   if(minimizer=="minuit"){
1489     res = new AliTrackResidualsChi2();
1490   }else if(minimizer=="minuitnorot"){
1491     res = new AliTrackResidualsChi2();
1492     res->FixParameter(3);
1493     res->FixParameter(4);
1494     res->FixParameter(5);
1495   }else if(minimizer=="fast"){
1496     res = new AliTrackResidualsFast();
1497   }else {
1498     printf("Trying to set a non existing minimizer! \n");
1499     return;
1500   }
1501
1502   res->SetMinNPoints(1);
1503   SetMinimizer(res);
1504   Bool_t draw = kTRUE;
1505
1506   if(misalignmentFile=="")printf("NO FAKE MISALIGNMENT\n");
1507   else {
1508     Bool_t misal=Misalign(misalignmentFile,"ITSAlignObjs");
1509     if(!misal)return;
1510   }
1511
1512   if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1513
1514   TArrayI volIDsFit(2200);
1515   Int_t iLayer,j=0;
1516   for (iLayer=(Int_t)AliGeomManager::kSPD1;iLayer<(Int_t)AliGeomManager::kTPC1;iLayer++){
1517     for (Int_t iModule=0;iModule<AliGeomManager::LayerSize(iLayer);iModule++){
1518       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iModule);
1519
1520       if((iLayer!=iLayerToAlign)&&(iLayer!=iLayerToExclude))volIDsFit.AddAt(volid,j);
1521
1522       j++;
1523     }
1524   }
1525
1526   Int_t size=AliGeomManager::LayerSize(iLayerToAlign);
1527   TArrayI volIDs(size);
1528
1529   j=0;
1530   for (Int_t iModule=0;iModule<AliGeomManager::LayerSize(iLayerToAlign);iModule++){
1531
1532     UShort_t volid = AliGeomManager::LayerToVolUID(iLayerToAlign,iModule);
1533     volIDs.AddAt(volid,j);
1534     j++;
1535   }
1536
1537     if(j==0){printf("j=0 \n");return;}
1538
1539     CalculateResiduals(&volIDs,&volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,1,draw);
1540
1541
1542     return;
1543
1544
1545 }
1546
1547 //______________________________________________________________________________
1548 void AliITSResidualsAnalysis::ExtractResiduals(Int_t layer,
1549                                                Int_t minEnt,
1550                                                TString filename) const
1551
1552 {
1553
1554   //
1555   // Function that saves the residuals into a Entuple
1556   //
1557
1558   TString title,strminEnt=" (Npts > ";
1559   histProperties_t histRPHIprop,histZprop;
1560
1561   // Labels for the plots
1562   strminEnt+=minEnt;
1563   strminEnt.Append(")");
1564
1565   // name of the output file
1566   title="resPlot_MA_layer";
1567   title+=layer;
1568   title.Append(".root");
1569
1570   // Load INfiles, OUTfiles and TTrees and labels them
1571   TFile *f1=TFile::Open(filename.Data());
1572   TFile &f2=*f1;
1573   TFile *outfile2=new TFile(title.Data(),"RECREATE");
1574   TFile &outfile=*outfile2;
1575   TTree *tRealign2=(TTree*)f2.Get("analysisTree"); // TTree with the DATA
1576   TTree &tRealign=*tRealign2;
1577
1578
1579   // Setting variables
1580   Int_t nEntries;
1581   Int_t *volid;
1582   Float_t z,phi;
1583   TH2F *hVolCorrBranch;
1584   TH1F *hResRPhi;
1585   TH1F *hResZ;
1586
1587   TString layer2="";
1588   layer2+=layer;
1589
1590
1591   TH2F *hEmpty=(TH2F*)f2.Get("fhEmpty");
1592   hEmpty->SetName("hEmpty");
1593
1594
1595   // Creates histos using the "hEmpty" template (binning!)
1596   TH2F *hSigmaMeanRPHI=new TH2F();
1597   TH2F *hSigmaRPHI=new TH2F();
1598   TH2F *hSigmaMeanZ=new TH2F();
1599   hEmpty->Copy(*hSigmaMeanRPHI);
1600   hSigmaMeanRPHI->SetName("hSigmaMeanRPHI");
1601   hSigmaMeanRPHI->GetZaxis()->SetRangeUser(0.,200);
1602   hEmpty->Copy(*hSigmaRPHI);
1603   hSigmaRPHI->SetName("hSigmaRPHI");
1604   hSigmaRPHI->GetZaxis()->SetRangeUser(0.,200);
1605   hEmpty->Copy(*hSigmaMeanZ);
1606   hSigmaMeanZ->SetName("hSigmaMeanZ");
1607   hSigmaMeanZ->GetZaxis()->SetRangeUser(0.,400);
1608
1609
1610   // Branching of the DATA TTree
1611   tRealign.SetBranchAddress("globalPhi",&phi);
1612   tRealign.SetBranchAddress("globalZ",&z);
1613   tRealign.SetBranchAddress("histZ",&hResZ);
1614   tRealign.SetBranchAddress("histRPHI",&hResRPhi);
1615   tRealign.SetBranchAddress("volID",&volid);
1616   tRealign.SetBranchAddress("histCorrVol",&hVolCorrBranch);
1617   tRealign.SetBranchAddress("histRPHIprop",&histRPHIprop);
1618   tRealign.SetBranchAddress("histZprop",&histZprop);
1619
1620   TNtuple *ntMonA = new TNtuple("ntMonA","Residuals","layer:volID:phi:z:nentries:meanFitRPHI:meanFitZ:RMS_RPHI");
1621   nEntries=tRealign.GetEntries();
1622   printf("entries: %d\n",nEntries);
1623   Float_t volidfill = 0;
1624
1625   for(Int_t j=0;j<nEntries;j++){ // LOOP OVER THE ENTRIES
1626
1627     printf(" Loading Event %d \n",j);
1628
1629     tRealign.GetEvent(j);
1630
1631     // Saving data in an entuple -> To be turned into a Tree
1632     ntMonA->Fill((Float_t)layer,
1633                 volidfill, // WRONG! To be corrected!
1634                 (Float_t)phi,
1635                 (Float_t)z,
1636                 10000*(Float_t)histRPHIprop.nentries,
1637                 10000*(Float_t)histRPHIprop.meanFit,
1638                 10000*(Float_t)histZprop.meanFit,
1639                 10000*(Float_t)histRPHIprop.rms);
1640
1641   } // END LOOP OVER ENTRIES
1642
1643
1644   //write histograms
1645   outfile.cd();//return to top directory
1646   ntMonA->Write();
1647   hEmpty->Write();
1648
1649   delete  tRealign2;
1650   f2.Close();
1651
1652   return;
1653
1654 }
1655
1656 //______________________________________________________________________________
1657 Int_t AliITSResidualsAnalysis::PlotResiduals(Int_t layer,TString filename) const
1658 {
1659   //
1660   // Function that plots the residual distributions
1661   //
1662   filename+=layer;
1663   filename.Append(".root");
1664   TFile *f1 = TFile::Open(filename.Data());
1665   if(!f1) return 1;
1666
1667   TH2F *hEmpty=(TH2F*)f1->Get("hEmpty");
1668   TNtuple *ntData=(TNtuple*)f1->Get("ntMonA");
1669   if(!ntData) return 2;
1670
1671
1672   TH2F *hMeanZ = new TH2F("hMeanZ","Hist for getting banged",40,-20,20,30,-15,15);
1673
1674
1675   Int_t nn=4;
1676   Float_t x[4],y[4],ex[4],ey[4],yZ[4],eyZ[4];
1677   x[0]=10.5;
1678   x[1]=3.5;
1679   x[2]=-3.5;
1680   x[3]=-10.5;
1681
1682   // Declaring Histos
1683   TH2F *hStatGlob = new TH2F();
1684   TH2F *hMeanGlob = new TH2F();
1685
1686   TH1F **hMeanPHI;
1687   TH1F **hMeanPHIz;
1688   TH1F *hGlobPhi = new TH1F("hGlobPhi","hGlobPhi",41,-(TMath::Pi())-(TMath::Pi()/40),(TMath::Pi())+(TMath::Pi()/40));
1689   //TH1F *hGlobPhi = new TH1F("hGlobPhi","hGlobPhi",40,-(TMath::Pi()),(TMath::Pi()));
1690
1691   hMeanPHI = new TH1F*[40]; //watch out!
1692   hMeanPHIz = new TH1F*[40];
1693
1694   TString title;
1695
1696   for(Int_t bPhi = 0; bPhi<40; bPhi++){
1697     title="hMeanPHI";
1698     title+=bPhi;
1699     hMeanPHI[bPhi]=new TH1F(title.Data(),title.Data(),300,-150,150);
1700     title="hMeanPHIz";
1701     title+=bPhi;
1702     hMeanPHIz[bPhi]=new TH1F(title.Data(),title.Data(),300,-150,150);
1703   }
1704
1705   // Setting the binning of the histos
1706   hEmpty->Copy(*hStatGlob);
1707   hStatGlob->SetName("hStatGlob");
1708   hEmpty->Copy(*hMeanGlob);
1709   hMeanGlob->SetName("hMeanGlob");
1710
1711   Int_t entries;
1712   Float_t volID,phi,z,nentries,meanFitRPHI,meanFitZ,rms;
1713   entries = (Int_t)ntData->GetEntries();
1714
1715   // Branching ...
1716   //ntData->SetBranchAddress("layer",&layernt);
1717   ntData->SetBranchAddress("volID",&volID);
1718   ntData->SetBranchAddress("phi",&phi);
1719   ntData->SetBranchAddress("z",&z);
1720   ntData->SetBranchAddress("nentries",&nentries);
1721   ntData->SetBranchAddress("meanFitRPHI",&meanFitRPHI);
1722   ntData->SetBranchAddress("meanFitZ",&meanFitZ);
1723   ntData->SetBranchAddress("RMS_RPHI",&rms);
1724
1725   Int_t nbytes = 0;
1726   Int_t bin,bin2,ban;
1727   Double_t c1,c2,c3,c4;
1728   Double_t m1,m2,m3,m4;
1729   Double_t n1,n2,n3,n4;
1730   c1=1e-10;
1731   c2=1e-10;
1732   c3=1e-10;
1733   c4=1e-10;
1734
1735   TH1F *hMeanFit1 = new TH1F("hMeanFit1","lol",1000,-500,500);
1736   TH1F *hMeanFit2 = new TH1F("hMeanFit2","lol",1000,-500,500);
1737   TH1F *hMeanFit3 = new TH1F("hMeanFit3","lol",1000,-500,500);
1738   TH1F *hMeanFit4 = new TH1F("hMeanFit4","lol",1000,-500,500);
1739
1740   TH1F *hMeanFitZ1 = new TH1F("hMeanFitZ1","lol",1000,-500,500);
1741   TH1F *hMeanFitZ2 = new TH1F("hMeanFitZ2","lol",1000,-500,500);
1742   TH1F *hMeanFitZ3 = new TH1F("hMeanFitZ3","lol",1000,-500,500);
1743   TH1F *hMeanFitZ4 = new TH1F("hMeanFitZ4","lol",1000,-500,500);
1744
1745   for(Int_t j=0;j<entries;j++){
1746
1747     nbytes += ntData->GetEvent(j);
1748
1749     // Set binning
1750     bin=hStatGlob->FindBin(phi,z);
1751     bin2=hMeanZ->FindBin(meanFitRPHI,z);
1752
1753     // Global Histograms
1754     hStatGlob->AddBinContent(bin,nentries);
1755     hMeanGlob->AddBinContent(bin,meanFitRPHI);
1756     hMeanZ->AddBinContent(bin2,1);
1757
1758     bin=hGlobPhi->FindBin(phi);
1759     bin2=hMeanPHI[bin-2]->FindBin(meanFitRPHI);
1760     hMeanPHI[bin-2]->AddBinContent(bin2,1);
1761     bin2=hMeanPHIz[bin-2]->FindBin(meanFitZ);
1762     hMeanPHIz[bin-2]->AddBinContent(bin2,1);
1763
1764
1765     if(z<12 && z>9) {
1766       c1+=nentries;
1767       m1+=(meanFitRPHI*nentries);
1768       n1+=rms*nentries;
1769       ban=hMeanFit1->FindBin(meanFitRPHI);
1770       //hMeanFit1->AddBinContent(ban,1);
1771       hMeanFit1->AddBinContent(ban,1);
1772       ban=hMeanFitZ1->FindBin(meanFitZ);
1773       hMeanFitZ1->AddBinContent(ban,1);
1774     } else if(z<5 && z>2){
1775       c2+=nentries;
1776       m2+=(meanFitRPHI*nentries);
1777       n2+=rms*nentries;
1778       ban=hMeanFit2->FindBin(meanFitRPHI);
1779       //hMeanFit2->AddBinContent(ban,1);
1780       hMeanFit2->AddBinContent(ban,1);
1781       ban=hMeanFitZ2->FindBin(meanFitZ);
1782       hMeanFitZ2->AddBinContent(ban,1);
1783     } else if(z<-2 && z>-5){
1784       c3+=nentries;
1785       m3+=(meanFitRPHI*nentries);
1786       n3+=rms*nentries;
1787       ban=hMeanFit3->FindBin(meanFitRPHI);
1788       //hMeanFit3->AddBinContent(ban,1);
1789       hMeanFit3->AddBinContent(ban,1);
1790       ban=hMeanFitZ3->FindBin(meanFitZ);
1791       hMeanFitZ3->AddBinContent(ban,1);
1792     } else if(z<-9 && z>-12){
1793       c4+=nentries;
1794       m4+=(meanFitRPHI*nentries);
1795       n4+=rms*nentries;
1796       ban=hMeanFit4->FindBin(meanFitRPHI);
1797       //hMeanFit4->AddBinContent(ban,1);
1798       hMeanFit4->AddBinContent(ban,1);
1799       ban=hMeanFitZ4->FindBin(meanFitZ);
1800       hMeanFitZ4->AddBinContent(ban,1);
1801     }
1802
1803   }
1804
1805   ex[0]=3;
1806   ex[1]=3;
1807   ex[2]=3;
1808   ex[3]=3;
1809
1810   y[0]=hMeanFit1->GetMean();
1811   y[1]=hMeanFit2->GetMean();
1812   y[2]=hMeanFit3->GetMean();
1813   y[3]=hMeanFit4->GetMean();
1814
1815   ey[0]=hMeanFit1->GetRMS();
1816   ey[1]=hMeanFit2->GetRMS();
1817   ey[2]=hMeanFit3->GetRMS();
1818   ey[3]=hMeanFit4->GetRMS();
1819
1820   yZ[0]=hMeanFitZ1->GetMean();
1821   yZ[1]=hMeanFitZ2->GetMean();
1822   yZ[2]=hMeanFitZ3->GetMean();
1823   yZ[3]=hMeanFitZ4->GetMean();
1824
1825   eyZ[0]=hMeanFitZ1->GetRMS();
1826   eyZ[1]=hMeanFitZ2->GetRMS();
1827   eyZ[2]=hMeanFitZ3->GetRMS();
1828   eyZ[3]=hMeanFitZ4->GetRMS();
1829
1830   TGraphErrors *gZres = new TGraphErrors(nn,x,y,ex,ey);
1831   TGraphErrors *gZresZ = new TGraphErrors(nn,x,yZ,ex,eyZ);
1832
1833   TCanvas *cc1 = new TCanvas("cc1","Title1",1);
1834   cc1->cd();
1835   hStatGlob->DrawCopy("LEGO2");
1836
1837   Double_t xx[40],yy[40],exx[40],eyy[40];
1838
1839   for(Int_t bp = 0; bp<40;bp++){
1840     xx[bp]=(bp*(TMath::Pi()/20))-TMath::Pi();
1841     if(TMath::Abs(hMeanPHI[bp]->GetMean())<1e-6) continue;
1842     yy[bp]=hMeanPHI[bp]->GetMean();
1843     eyy[bp]=hMeanPHI[bp]->GetRMS();
1844     exx[bp]=(TMath::Pi())/41;
1845   }
1846   TGraphErrors *gPHIres = new TGraphErrors(40,xx,yy,exx,eyy);
1847   //gPHIres->Fit("pol1","","same",-3,3);
1848
1849   Double_t xxz[40],yyz[40],exxz[40],eyyz[40];
1850
1851   for(Int_t bp = 0; bp<40;bp++){
1852     xxz[bp]=(bp*(TMath::Pi()/20))-TMath::Pi();
1853     if(TMath::Abs(hMeanPHIz[bp]->GetMean())<1e-6) continue;
1854     yyz[bp]=hMeanPHIz[bp]->GetMean();
1855     eyyz[bp]=hMeanPHIz[bp]->GetRMS();
1856     exxz[bp]=(TMath::Pi())/41;
1857   }
1858   TGraphErrors *gPHIresZ = new TGraphErrors(40,xxz,yyz,exxz,eyyz);
1859
1860   TCanvas *cc4 = new TCanvas("cc4","meanRes VS Z",1);
1861   cc4->Divide(1,2);
1862   cc4->cd(1);
1863   gZres->Draw("AP");
1864   cc4->cd(2);
1865   gZresZ->Draw("AP");
1866
1867   TCanvas *cc3 = new TCanvas("cc3","Title3",1);
1868   cc3->Divide(2,2);
1869   cc3->cd(1);
1870   hMeanFitZ1->DrawCopy();
1871   cc3->cd(2);
1872   hMeanFitZ2->DrawCopy();
1873   cc3->cd(3);
1874   hMeanFitZ3->DrawCopy();
1875   cc3->cd(4);
1876   hMeanFitZ4->DrawCopy();
1877
1878   TCanvas *cc6 = new TCanvas("cc6","meanRes(RPHI) VS PHI",1);
1879   cc6->Divide(1,2);
1880   cc6->cd(1);
1881   gPHIres->Draw("AP");
1882
1883   cc6->cd(2);
1884   gPHIresZ->Draw("AP");
1885
1886   TCanvas *cc7 = new TCanvas("cc7","Title7",1);
1887   cc7->Divide(2,2);
1888   cc7->cd(1);
1889   hMeanPHI[1]->DrawCopy();
1890   cc7->cd(2);
1891   hMeanPHI[2]->DrawCopy();
1892   cc7->cd(3);
1893   hMeanPHI[29]->DrawCopy();
1894   cc7->cd(4);
1895   hMeanPHI[30]->DrawCopy();
1896
1897
1898
1899   f1->Close();
1900
1901   return 0;
1902 }
1903