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