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