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