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