1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //-----------------------------------------------------------------
25 #include <Riostream.h>
30 #include <TClonesArray.h>
37 #include <TGraphErrors.h>
38 #include "TGeoMatrix.h"
39 #include "TGeoManager.h"
40 #include "TGeoPhysicalNode.h"
41 #include "TMatrixDSymEigen.h"
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"
56 #include "AliITSResidualsAnalysis.h"
59 // Structure of the RealignmentAnalysisLayer.C
68 ClassImp(AliITSResidualsAnalysis)
70 //____________________________________________________________________________
71 AliITSResidualsAnalysis::AliITSResidualsAnalysis():
88 fTrackDirLambdaAll(0),
89 fTrackDirLambda2All(0),
107 fRealignObjFileIsOpen(kFALSE),
109 fAliTrackPoints("AliTrackPoints.root"),
110 fGeom("geometry.root")
121 //____________________________________________________________________________
122 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TString aliTrackPoints,
124 AliAlignmentTracks(),
140 fTrackDirLambdaAll(0),
141 fTrackDirLambda2All(0),
142 fTrackDirAlphaAll(0),
159 fRealignObjFileIsOpen(kFALSE),
161 fAliTrackPoints(aliTrackPoints),
165 // Constructor (alitrackpoints)
171 //____________________________________________________________________________
172 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TArrayI *volIDs):
173 AliAlignmentTracks(),
189 fTrackDirLambdaAll(0),
190 fTrackDirLambda2All(0),
191 fTrackDirAlphaAll(0),
208 fRealignObjFileIsOpen(kFALSE),
210 fAliTrackPoints("AliTrackPoints.root"),
211 fGeom("geometry.root")
215 // Original Constructor
218 InitHistograms(volIDs);
222 //____________________________________________________________________________
223 AliITSResidualsAnalysis::AliITSResidualsAnalysis(TArrayI *volIDs,AliTrackPointArray **tracksClustArray,AliTrackPointArray **tracksFitPointsArray):
224 AliAlignmentTracks(),
240 fTrackDirLambdaAll(0),
241 fTrackDirLambda2All(0),
242 fTrackDirAlphaAll(0),
259 fRealignObjFileIsOpen(kFALSE),
261 fAliTrackPoints("AliTrackPoints.root"),
262 fGeom("geometry.root")
266 // Detailed Constructor (deprecated)
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);
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 ");
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 ");
319 fpTrackVolIDs=new TArrayI(fnHist);
320 fVolUsed=new TArrayI*[fnHist];
321 fVolVolids=new TArrayI*[fnHist];
322 fLastVolVolid=new Int_t[fnHist];
324 for (Int_t nhist=0;nhist<fnHist;nhist++){
325 fpTrackVolIDs->AddAt(volIDs->At(nhist),nhist);
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());
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());
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());
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());
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());
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 ");
364 fVolVolids[nhist]=new TArrayI(1000);
365 fVolUsed[nhist]=new TArrayI(1000);
366 fLastVolVolid[nhist]=0;
367 FillResHists(tracksClustArray[nhist],tracksFitPointsArray[nhist]);
372 SetFileNameTrackPoints(""); // Filename with the AliTrackPoints
373 SetFileNameGeometry(""); // Filename with the Geometry
378 //____________________________________________________________________________
379 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const AliITSResidualsAnalysis& /* obj */):
380 AliAlignmentTracks(),
396 fTrackDirLambdaAll(0),
397 fTrackDirLambda2All(0),
398 fTrackDirAlphaAll(0),
415 fRealignObjFileIsOpen(kFALSE),
417 fAliTrackPoints("AliTrackPoints.root"),
418 fGeom("geometry.root")
421 // copy constructor. This is not allowed.
423 AliFatal("Copy constructor not allowed\n");
427 //____________________________________________________________________________
428 AliITSResidualsAnalysis& AliITSResidualsAnalysis::operator = (const AliITSResidualsAnalysis& /* obj */) {
429 // assignment operator. This is not allowed
430 AliFatal("Assignment operator not allowed\n");
434 //____________________________________________________________________________
435 AliITSResidualsAnalysis::~AliITSResidualsAnalysis()
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;
472 SetFileNameTrackPoints("");
473 SetFileNameGeometry("");
477 //____________________________________________________________________________
478 void AliITSResidualsAnalysis::InitHistograms(const TArrayI *volIDs)
481 // Method that sets and creates the required hisstrograms
482 // with the correct binning (it dos not fill them)
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_";
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_";
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];
508 fTrackDirPhi=new TH1F*[fnHist];
509 fTrackDirLambda=new TH1F*[fnHist];
510 fTrackDirLambda2=new TH1F*[fnHist];
511 fTrackDirAlpha=new TH1F*[fnHist];
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);
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.);
523 fTrackDir=new TH2F*[fnHist];
525 Float_t **binningZPhi=CheckSingleLayer(volIDs);
526 fvolidsToBin=new Int_t*[fnPhi*fnZ];
528 Float_t *binningphi=binningZPhi[0];
529 Float_t *binningz=binningZPhi[1];
530 Bool_t binning=SetBinning(volIDs,binningphi,binningz);
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 ");
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 ");
562 fpTrackVolIDs=new TArrayI(fnHist);
563 fVolUsed=new TArrayI*[fnHist];
564 fVolVolids=new TArrayI*[fnHist];
565 fLastVolVolid=new Int_t[fnHist];
567 for (Int_t nhist=0;nhist<fnHist;nhist++){
568 fpTrackVolIDs->AddAt(volIDs->At(nhist),nhist);
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());
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());
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());
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());
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());
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());
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());
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());
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());
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());
630 aux+=volIDs->At(nhist);
631 fhistCorrVol[nhist]=new TH2F("histname","histname",50,-3.2,3.2,100,-80,80);
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;
648 //____________________________________________________________________________
649 void AliITSResidualsAnalysis::ListVolUsed(TTree *pointsTree,TArrayI ***arrayIndex,Int_t **lastIndex)
652 // This is copied from AliAlignmentClass::LoadPoints() method
655 Int_t volIDalignable,volIDpoint,iModule;
657 AliTrackPointArray* array = 0;
658 pointsTree->SetBranchAddress("SP", &array);
661 for(Int_t ivol=0;ivol<fnHist;ivol++){
663 volIDalignable=fpTrackVolIDs->At(ivol);
664 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volIDalignable,iModule);
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++) {
672 Int_t entry = (*index)[iArrayId];
674 pointsTree->GetEvent(entry);
676 AliWarning("Wrong space point array index!");
680 // Get the space-point array
681 Int_t modnum,nPoints = array->GetNPoints();
683 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
684 array->GetPoint(p,iPoint);
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));
694 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
696 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
697 layer,modnum,AliGeomManager::LayerSize(layer)));
700 if (layer > AliGeomManager::kSSD2) continue; // ITS only
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);
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){
719 size = fVolUsed[ivol]->GetSize();
720 // If needed allocate new size
721 if (lastused>= size){
722 fVolUsed[ivol]->Set(size + 1000);
724 fVolUsed[ivol]->AddAt(volIDpoint,lastused);
728 FillVolumeCorrelationHists(ivol,volIDalignable,volIDpoint,usedVol);
736 //____________________________________________________________________________
737 void AliITSResidualsAnalysis::FillVolumeCorrelationHists(Int_t ivol,Int_t volIDalignable,Int_t volIDpoint,Bool_t usedVol) const
740 // Fill the histograms with the correlations between volumes
743 if(!gGeoManager)AliGeomManager::LoadGeometry(GetFileNameGeometry());
744 Double_t *transGlobal,radius,phi;
745 const char *symname,*volpath;
747 TGeoPhysicalNode *pn;
748 TGeoHMatrix *globMatrix;
750 symname = AliGeomManager::SymName(volIDalignable);
751 pne = gGeoManager->GetAlignableEntry(symname);
752 volpath=pne->GetTitle();
753 pn=gGeoManager->MakePhysicalNode(volpath);
754 globMatrix=pn->GetMatrix();
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]);
762 symname = AliGeomManager::SymName(volIDpoint);
763 pne = gGeoManager->GetAlignableEntry(symname);
764 volpath=pne->GetTitle();
765 pn=gGeoManager->MakePhysicalNode(volpath);
766 globMatrix=pn->GetMatrix();
768 transGlobal=globMatrix->GetTranslation();
769 radius=TMath::Sqrt(transGlobal[0]*transGlobal[0]+transGlobal[1]*transGlobal[1]);
770 phi=TMath::ATan2(transGlobal[1],transGlobal[0]);
772 fhistCorrVol[ivol]->Fill(phi,transGlobal[2]);
778 //____________________________________________________________________________
779 void AliITSResidualsAnalysis::FillResHists(AliTrackPointArray *points,AliTrackPointArray *pTrack) const
782 // Method that fills the histograms with the residuals
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;
792 for(Int_t ipoint=0;ipoint<points->GetNPoints();ipoint++){
793 points->GetPoint(p,ipoint);
794 volIDpoint=(Int_t)p.GetVolumeID();
797 pTrack->GetPoint(pTr,ipoint);
798 GetTrackDirClusterCov(&pTr,phi,lambda,lambda2,alpha,xovery,zovery);
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]);
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]));
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);
824 fTrackDirPhi[ivolIDs]->Fill(phi);
825 fTrackDirLambda[ivolIDs]->Fill(lambda);
826 fTrackDirLambda2[ivolIDs]->Fill(lambda2);
827 fTrackDirAlpha[ivolIDs]->Fill(alpha);
829 fTrackDirPhiAll->Fill(phi);
830 fTrackDirLambdaAll->Fill(lambda);
831 fTrackDirLambda2All->Fill(lambda2);
832 fTrackDirAlphaAll->Fill(alpha);
834 fTrackDirAll->Fill(lambda,alpha);
835 fTrackDir2All->Fill(lambda2,phi);
836 fTrackDirXZAll->Fill(xovery,zovery);
837 fTrackDir[ivolIDs]->Fill(lambda,alpha);
839 fPullHistRPHI[ivolIDs]->Fill(pullrphi);
840 fPullHistZ[ivolIDs]->Fill(pullz);
844 Float_t globalPhi,globalZ;
845 if(kTRUE||(fvolidsToBin[ivolIDs][0]!=volIDpoint)){
846 binphi=GetBinPhiZ((Int_t)volIDpoint,&binz);
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];
852 globalPhi=fCoordToBinTable[binphi][binz][0];
853 globalZ=fCoordToBinTable[binphi][binz][1];
855 fVolNTracks->Fill(globalPhi,globalZ);
857 else fVolNTracks->Fill(TMath::ATan2(xyz[1],xyz[0]),xyz[2]);
864 //____________________________________________________________________________
865 Bool_t AliITSResidualsAnalysis::AnalyzeHists(Int_t minNpoints) const
868 // Saves the histograms into a tree and saves the tree into a file
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;
878 TH1F *histRPHI,*histZ,*histGlob,*histPullRPHI,*histPullZ,*hTrackDirPhi,*hTrackDirLambda,*hTrackDirLambda2,*hTrackDirAlpha;
880 TH2F *histCorrVol,*hTrackDir;
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();
892 histCorrVol=new TH2F();
893 Float_t globalPhi,globalZ;
895 Int_t nHistAnalyzed=0,entries;
896 analysisTree->Branch("volID",&volID,"volID/I");
898 analysisTree->Branch("globalPhi",&globalPhi,"globalPhi/F");
899 analysisTree->Branch("globalZ",&globalZ,"globalZ/F");
901 analysisTree->Branch("histRPHI","TH1F",&histRPHI,128000,0);
902 analysisTree->Branch("histPullRPHI","TH1F",&histPullRPHI,128000,0);
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);
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");
917 analysisTree->Branch("histCorrVol","TH2F",&histCorrVol,128000,0);
920 for(Int_t j=0;j<fnHist;j++){
921 volID=fpTrackVolIDs->At(j);
922 if((entries=(fResHistGlob[j]->GetEntries())>=minNpoints)||fsingleLayer){
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);
937 histPullZ=fPullHistZ[j];
938 histZprop.nentries=(Int_t)fResHistZ[j]->GetEntries();
939 rms=fResHistZ[j]->GetRMS();
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);
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);
958 hTrackDirPhi=fTrackDirPhi[j];
959 hTrackDirLambda=fTrackDirLambda[j];
960 hTrackDirLambda2=fTrackDirLambda2[j];
961 hTrackDirAlpha=fTrackDirAlpha[j];
962 hTrackDir=fTrackDir[j];
966 if (fvolidsToBin[j][0]!=volID)binphi=GetBinPhiZ((Int_t)volID,&binz);
968 binphi=fvolidsToBin[j][1];
969 binz=fvolidsToBin[j][2];
971 globalPhi=fCoordToBinTable[binphi][binz][0];
972 globalZ=fCoordToBinTable[binphi][binz][1];
975 histCorrVol=fhistCorrVol[j];
976 fSigmaVolZ->SetBinContent(binphi+1,binz+1,histRPHIprop.sigmaFit);//+1 is for underflow
978 analysisTree->Fill();
984 analysisTree->Print();
985 fVolNTracks->Write();
989 fhistVolUsed->Draw();
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();
1006 delete analysisTree;
1012 //____________________________________________________________________________
1013 void AliITSResidualsAnalysis::DrawHists() const
1016 // Draws the histograms of the residuals and of the number of tracks
1020 for(Int_t canv=0;canv<fnHist;canv++){
1021 cname="canv Residuals";
1023 TCanvas *c=new TCanvas(cname.Data(),cname.Data(),700,700);
1026 fVolResHistRPHI[canv]->Draw();
1028 fResHistZ[canv]->Draw();
1030 fResHistGlob[canv]->Draw();
1032 cname="canv NVolTracks";
1034 TCanvas *c2=new TCanvas(cname.Data(),cname.Data(),700,700);
1036 fVolNTracks->Draw();
1041 //____________________________________________________________________________
1042 Float_t** AliITSResidualsAnalysis::CheckSingleLayer(const TArrayI *volids)
1045 // Checks if volumes array is a single (ITS) layer or not
1048 Float_t **binningzphi=new Float_t*[2];
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;
1060 //Bool_t used=kFALSE;
1062 case AliGeomManager::kSPD1:{
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];
1068 for(Int_t j=0;j<fnPhi;j++){
1069 fCoordToBinTable[j]=new Double_t*[kZSPD1];
1072 case AliGeomManager::kSPD2:{
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];
1078 for(Int_t j=0;j<fnPhi;j++){
1079 fCoordToBinTable[j]=new Double_t*[kZSPD2];
1082 }; break; case AliGeomManager::kSDD1:{
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];
1088 for(Int_t j=0;j<fnPhi;j++){
1089 fCoordToBinTable[j]=new Double_t*[kZSDD1];
1091 }; break; case AliGeomManager::kSDD2:{
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];
1097 for(Int_t j=0;j<fnPhi;j++){
1098 fCoordToBinTable[j]=new Double_t*[kZSDD2];
1100 }; break; case AliGeomManager::kSSD1:{
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];
1106 for(Int_t j=0;j<fnPhi;j++){
1107 fCoordToBinTable[j]=new Double_t*[kZSSD1];
1109 }; break; case AliGeomManager::kSSD2:{
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];
1115 for(Int_t j=0;j<fnPhi;j++){
1116 fCoordToBinTable[j]=new Double_t*[kZSSD2];
1121 printf("Wrong Layer Label! \n");
1122 fsingleLayer=kFALSE;
1130 //____________________________________________________________________________
1131 Bool_t AliITSResidualsAnalysis::SetBinning(const TArrayI *volids,Float_t *phiBin,Float_t *zBin)
1134 // Sets the correct binning
1137 if(!fsingleLayer)return kFALSE;
1138 const char *volpath,*symname;
1140 Int_t *orderArrayPhi,*orderArrayZ;
1142 Double_t *phiArray,*zArray,*transGlobal,*phiArrayOrdered,*zArrayOrdered;
1143 Double_t lastPhimin=-10;
1144 Double_t lastZmin=-99;
1147 TGeoPhysicalNode *pn;
1148 TGeoHMatrix *globMatrix;
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++){
1163 for(Int_t k=0;k<fnPhi;k++){
1166 orderPhiZ[k]=new Int_t*[fnZ];
1170 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volids->At(0),iModule);
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();
1184 for(Int_t j=0;j<lastPhi;j++){
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;
1193 phiArray[lastPhi]=TMath::ATan2(transGlobal[1],transGlobal[0]);
1194 fvolidsToBin[iModule][1]=lastPhi;
1195 if(phiArray[lastPhi]<lastPhimin)lastPhimin=phiArray[lastPhi];
1198 printf("Wrong Phi! \n");
1202 for(Int_t j=0;j<lastZ;j++){
1204 if(TMath::Abs(zArray[j]-transGlobal[2])<0.1){
1205 fvolidsToBin[iModule][2]=j;
1211 fvolidsToBin[iModule][2]=lastZ;
1212 zArray[lastZ]=transGlobal[2];
1213 if(zArray[lastZ]<lastZmin)lastZmin=zArray[lastZ];
1216 printf("Wrong Z! \n");
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]++;
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];
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]++;
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];
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];
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]);
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];
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]);
1288 phiBin[0]=(Float_t)phiArrayOrdered[0]-(phiArrayOrdered[1]-phiArrayOrdered[0])/2.;
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]);
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]);
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]);
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]);
1309 //____________________________________________________________________________
1310 Int_t AliITSResidualsAnalysis::GetBinPhiZ(const Int_t volID,Int_t *binz) const
1313 // Returns the correct Phi-Z bin
1317 printf("No Single Layer reAlignment! \n");
1321 for(Int_t j=0;j<fnPhi*fnZ;j++){
1323 printf("Wrong volume UID introduced! fnHist: %d volID: %d iter: %d \n",fnHist,volID,j);
1326 if(fvolidsToBin[j][0]==volID){
1328 *binz=fvolidsToBin[j][2];
1329 return fvolidsToBin[j][1];
1337 //____________________________________________________________________________
1338 TArrayI* AliITSResidualsAnalysis::GetSingleLayerVolids(Int_t layer) const
1341 // Translates the layer number into a Volumes Array
1344 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1346 if(layer<1 || layer>6){
1347 printf("WRONG LAYER SET! \n");
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);
1363 //____________________________________________________________________________
1364 void AliITSResidualsAnalysis::GetTrackDirClusterCov(AliTrackPoint *point,Double_t &phi,Double_t &lambda,Double_t &lambda2,Double_t &alpha,Double_t &xovery,Double_t &zovery) const
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];
1379 Double_t determinant=cov.Determinant();
1380 if(determinant!=0.){
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];
1399 mainvect[0]=-1.*mainvect[0];
1400 mainvect[1]=-1.*mainvect[1];
1401 mainvect[2]=-1.*mainvect[2];
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();
1408 else printf("determinant =0!, skip this point \n");
1413 //____________________________________________________________________________
1414 void AliITSResidualsAnalysis::CalculateResiduals(const TArrayI *volids,
1415 const TArrayI *volidsfit,
1416 AliGeomManager::ELayerID layerRangeMin,
1417 AliGeomManager::ELayerID layerRangeMax,
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
1429 Int_t nVolIds = volids->GetSize();
1431 AliError("Volume IDs array is empty!");
1435 // Load only the tracks with at least one
1436 // space point in the set of volume (volids)
1438 //AliAlignmentTracks::SetPointsFilename(GetFileNameTrackPoints());
1439 AliAlignmentTracks::BuildIndex();
1441 cout<<" Index Built!"<<endl;
1443 if(draw) ListVolUsed(fPointsTree,fArrayIndex,fLastIndex);
1445 AliTrackPointArray **points;
1447 // Start the iterations
1448 while (iterations > 0) {
1449 Int_t nArrays = LoadPoints(volids, points);
1450 if (nArrays == 0) return;
1452 AliTrackResiduals *minimizer = CreateMinimizer();
1453 minimizer->SetNTracks(nArrays);
1454 minimizer->InitAlignObj();
1455 AliTrackFitter *fitter = CreateFitter();
1457 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
1458 if (!points[iArray]) continue;
1461 fitter->SetTrackPointArray(points[iArray],kFALSE);
1462 if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
1463 AliTrackPointArray *pVolId,*pTrack;
1466 fitter->GetTrackResiduals(pVolId,pTrack);
1467 if(draw) FillResHists(pVolId,pTrack); // WARNING!
1469 minimizer->AddTrackPointArrays(pVolId,pTrack);
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);
1479 minimizer->Minimize();
1481 // Update the alignment object(s)
1484 if(fRealignObjFileIsOpen)last=fClonesArray->GetLast();
1487 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
1488 UShort_t volid = (*volids)[iVolId];
1490 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
1491 AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
1492 *alignObj *= *minimizer->GetAlignObj();
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);
1504 UnloadPoints(nArrays, points);
1509 if(draw && iterations==0) AnalyzeHists(30);
1518 //______________________________________________________________________________
1519 void AliITSResidualsAnalysis::ProcessPoints(TString minimizer,
1521 AliGeomManager::ELayerID iLayerToAlign,
1522 AliGeomManager::ELayerID iLayerToExclude,
1523 TString misalignmentFile)
1526 // This function process the AliTrackPoints (into residuals)
1529 SetPointsFilename(GetFileNameTrackPoints());
1530 AliTrackFitter *fitter;
1533 fitter = new AliTrackFitterKalman();
1534 }else fitter = new AliTrackFitterRieman();
1536 fitter->SetMinNPoints(4);
1537 SetTrackFitter(fitter);
1540 AliTrackResiduals *res;
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();
1552 printf("Trying to set a non existing minimizer! \n");
1556 res->SetMinNPoints(1);
1558 Bool_t draw = kTRUE;
1560 if(misalignmentFile=="")printf("NO FAKE MISALIGNMENT\n");
1562 Bool_t misal=Misalign(misalignmentFile,"ITSAlignObjs");
1566 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1568 TArrayI volIDsFit(2200);
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);
1574 if((iLayer!=iLayerToAlign)&&(iLayer!=iLayerToExclude))volIDsFit.AddAt(volid,j);
1580 Int_t size=AliGeomManager::LayerSize(iLayerToAlign);
1581 TArrayI volIDs(size);
1584 for (Int_t iModule=0;iModule<AliGeomManager::LayerSize(iLayerToAlign);iModule++){
1586 UShort_t volid = AliGeomManager::LayerToVolUID(iLayerToAlign,iModule);
1587 volIDs.AddAt(volid,j);
1591 if(j==0){printf("j=0 \n");return;}
1593 CalculateResiduals(&volIDs,&volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,1,draw);
1601 //______________________________________________________________________________
1602 void AliITSResidualsAnalysis::ExtractResiduals(Int_t layer,
1604 TString filename) const
1609 // Function that saves the residuals into a Entuple
1612 TString title,strminEnt=" (Npts > ";
1613 histProperties_t histRPHIprop,histZprop;
1615 // Labels for the plots
1617 strminEnt.Append(")");
1619 // name of the output file
1620 title="resPlot_MA_layer";
1622 title.Append(".root");
1624 // Load INfiles, OUTfiles and TTrees and labels them
1625 TFile *f1=TFile::Open(filename.Data());
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;
1633 // Setting variables
1637 TH2F *hVolCorrBranch;
1645 TH2F *hEmpty=(TH2F*)f2.Get("fhEmpty");
1646 hEmpty->SetName("hEmpty");
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);
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);
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;
1679 for(Int_t j=0;j<nEntries;j++){ // LOOP OVER THE ENTRIES
1681 printf(" Loading Event %d \n",j);
1683 tRealign.GetEvent(j);
1685 // Saving data in an entuple -> To be turned into a Tree
1686 ntMonA->Fill((Float_t)layer,
1687 volidfill, // WRONG! To be corrected!
1690 10000*(Float_t)histRPHIprop.nentries,
1691 10000*(Float_t)histRPHIprop.meanFit,
1692 10000*(Float_t)histZprop.meanFit,
1693 10000*(Float_t)histRPHIprop.rms);
1695 } // END LOOP OVER ENTRIES
1699 outfile.cd();//return to top directory
1710 //______________________________________________________________________________
1711 Int_t AliITSResidualsAnalysis::PlotResiduals(Int_t layer,TString filename) const
1714 // Function that plots the residual distributions
1717 filename.Append(".root");
1718 TFile *f1 = TFile::Open(filename.Data());
1721 TH2F *hEmpty=(TH2F*)f1->Get("hEmpty");
1722 TNtuple *ntData=(TNtuple*)f1->Get("ntMonA");
1723 if(!ntData) return 2;
1726 TH2F *hMeanZ = new TH2F("hMeanZ","Hist for getting banged",40,-20,20,30,-15,15);
1730 Float_t x[4],y[4],ex[4],ey[4],yZ[4],eyZ[4];
1737 TH2F *hStatGlob = new TH2F();
1738 TH2F *hMeanGlob = new TH2F();
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()));
1745 hMeanPHI = new TH1F*[40]; //watch out!
1746 hMeanPHIz = new TH1F*[40];
1750 for(Int_t bPhi = 0; bPhi<40; bPhi++){
1753 hMeanPHI[bPhi]=new TH1F(title.Data(),title.Data(),300,-150,150);
1756 hMeanPHIz[bPhi]=new TH1F(title.Data(),title.Data(),300,-150,150);
1759 // Setting the binning of the histos
1760 hEmpty->Copy(*hStatGlob);
1761 hStatGlob->SetName("hStatGlob");
1762 hEmpty->Copy(*hMeanGlob);
1763 hMeanGlob->SetName("hMeanGlob");
1766 Float_t volID,phi,z,nentries,meanFitRPHI,meanFitZ,rms;
1767 entries = (Int_t)ntData->GetEntries();
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);
1781 Double_t c1,c2,c3,c4;
1782 Double_t m1,m2,m3,m4;
1783 Double_t n1,n2,n3,n4;
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);
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);
1799 for(Int_t j=0;j<entries;j++){
1801 nbytes += ntData->GetEvent(j);
1804 bin=hStatGlob->FindBin(phi,z);
1805 bin2=hMeanZ->FindBin(meanFitRPHI,z);
1807 // Global Histograms
1808 hStatGlob->AddBinContent(bin,nentries);
1809 hMeanGlob->AddBinContent(bin,meanFitRPHI);
1810 hMeanZ->AddBinContent(bin2,1);
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);
1821 m1+=(meanFitRPHI*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){
1830 m2+=(meanFitRPHI*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){
1839 m3+=(meanFitRPHI*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){
1848 m4+=(meanFitRPHI*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);
1864 y[0]=hMeanFit1->GetMean();
1865 y[1]=hMeanFit2->GetMean();
1866 y[2]=hMeanFit3->GetMean();
1867 y[3]=hMeanFit4->GetMean();
1869 ey[0]=hMeanFit1->GetRMS();
1870 ey[1]=hMeanFit2->GetRMS();
1871 ey[2]=hMeanFit3->GetRMS();
1872 ey[3]=hMeanFit4->GetRMS();
1874 yZ[0]=hMeanFitZ1->GetMean();
1875 yZ[1]=hMeanFitZ2->GetMean();
1876 yZ[2]=hMeanFitZ3->GetMean();
1877 yZ[3]=hMeanFitZ4->GetMean();
1879 eyZ[0]=hMeanFitZ1->GetRMS();
1880 eyZ[1]=hMeanFitZ2->GetRMS();
1881 eyZ[2]=hMeanFitZ3->GetRMS();
1882 eyZ[3]=hMeanFitZ4->GetRMS();
1884 TGraphErrors *gZres = new TGraphErrors(nn,x,y,ex,ey);
1885 TGraphErrors *gZresZ = new TGraphErrors(nn,x,yZ,ex,eyZ);
1887 TCanvas *cc1 = new TCanvas("cc1","Title1",1);
1889 hStatGlob->DrawCopy("LEGO2");
1891 Double_t xx[40],yy[40],exx[40],eyy[40];
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;
1900 TGraphErrors *gPHIres = new TGraphErrors(40,xx,yy,exx,eyy);
1901 //gPHIres->Fit("pol1","","same",-3,3);
1903 Double_t xxz[40],yyz[40],exxz[40],eyyz[40];
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;
1912 TGraphErrors *gPHIresZ = new TGraphErrors(40,xxz,yyz,exxz,eyyz);
1914 TCanvas *cc4 = new TCanvas("cc4","meanRes VS Z",1);
1921 TCanvas *cc3 = new TCanvas("cc3","Title3",1);
1924 hMeanFitZ1->DrawCopy();
1926 hMeanFitZ2->DrawCopy();
1928 hMeanFitZ3->DrawCopy();
1930 hMeanFitZ4->DrawCopy();
1932 TCanvas *cc6 = new TCanvas("cc6","meanRes(RPHI) VS PHI",1);
1935 gPHIres->Draw("AP");
1938 gPHIresZ->Draw("AP");
1940 TCanvas *cc7 = new TCanvas("cc7","Title7",1);
1943 hMeanPHI[1]->DrawCopy();
1945 hMeanPHI[2]->DrawCopy();
1947 hMeanPHI[29]->DrawCopy();
1949 hMeanPHI[30]->DrawCopy();