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