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