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>
29 #include <TClonesArray.h>
36 #include <TGraphErrors.h>
37 #include "TGeoMatrix.h"
38 #include "TGeoManager.h"
39 #include "TGeoPhysicalNode.h"
40 #include "TMatrixDSym.h"
41 #include "TMatrixDSymEigen.h"
45 #include "AliAlignmentTracks.h"
46 #include "AliTrackPointArray.h"
47 #include "AliAlignObjParams.h"
48 #include "AliTrackResiduals.h"
49 #include "AliTrackFitter.h"
50 #include "AliTrackFitterKalman.h"
51 #include "AliTrackFitterRieman.h"
52 #include "AliTrackResiduals.h"
53 #include "AliTrackResidualsChi2.h"
54 #include "AliTrackResidualsFast.h"
56 #include "AliITSgeomTGeo.h"
58 #include "AliITSResidualsAnalysis.h"
60 ClassImp(AliITSResidualsAnalysis)
62 //____________________________________________________________________________
63 AliITSResidualsAnalysis::AliITSResidualsAnalysis():
84 fTrackDirLambdaAll(0),
85 fTrackDirLambda2All(0),
103 fRealignObjFileIsOpen(kFALSE),
105 fAliTrackPoints("AliTrackPoints.root"),
106 fGeom("geometry.root"),
117 //____________________________________________________________________________
118 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TString aliTrackPoints,
120 AliAlignmentTracks(),
140 fTrackDirLambdaAll(0),
141 fTrackDirLambda2All(0),
142 fTrackDirAlphaAll(0),
159 fRealignObjFileIsOpen(kFALSE),
161 fAliTrackPoints(aliTrackPoints),
166 // Standard Constructor (alitrackpoints)
172 //____________________________________________________________________________
173 AliITSResidualsAnalysis::AliITSResidualsAnalysis(const TArrayI *volIDs):
174 AliAlignmentTracks(),
194 fTrackDirLambdaAll(0),
195 fTrackDirLambda2All(0),
196 fTrackDirAlphaAll(0),
213 fRealignObjFileIsOpen(kFALSE),
215 fAliTrackPoints("AliTrackPoints.root"),
216 fGeom("geometry.root"),
220 // Original Constructor
222 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
223 InitHistograms(volIDs);
227 //____________________________________________________________________________
228 AliITSResidualsAnalysis::~AliITSResidualsAnalysis()
234 if(fvolidsToBin) delete[] fvolidsToBin;
235 if(fLastVolVolid) delete[] fLastVolVolid;
236 if(fCoordToBinTable) delete[] fCoordToBinTable;
237 if(fHistCoordGlobY) delete[] fHistCoordGlobY;
238 if(fVolResHistRPHI) delete fVolResHistRPHI;
239 if(fResHistZ) delete fResHistZ;
241 for(Int_t i=0; i<fnHist; i++) delete fResHistX[i];
244 if(fResHistXLocsddL){
245 for(Int_t i=0; i<fnHist; i++) delete fResHistXLocsddL[i];
246 delete [] fResHistXLocsddL;
248 if(fResHistXLocsddR){
249 for(Int_t i=0; i<fnHist; i++) delete fResHistXLocsddR[i];
250 delete [] fResHistXLocsddR;
252 if(fPullHistRPHI) delete fPullHistRPHI;
253 if(fPullHistZ) delete fPullHistZ;
254 if(fTrackDirPhi) delete fTrackDirPhi;
255 if(fTrackDirLambda) delete fTrackDirLambda;
256 if(fTrackDirLambda2) delete fTrackDirLambda2;
257 if(fTrackDirAlpha) delete fTrackDirAlpha;
258 if(fTrackDirPhiAll) delete fTrackDirPhiAll;
259 if(fTrackDirLambdaAll) delete fTrackDirLambdaAll;
260 if(fTrackDirLambda2All) delete fTrackDirLambda2All;
261 if(fTrackDirAlphaAll) delete fTrackDirAlphaAll;
262 if(fTrackDir) delete fTrackDir;
263 if(fTrackDirAll) delete fTrackDirAll;
264 if(fTrackDir2All) delete fTrackDir2All;
265 if(fTrackDirXZAll) delete fTrackDirXZAll;
266 if(fResHistGlob) delete fResHistGlob;
267 if(fhistCorrVol) delete fhistCorrVol;
268 if(fVolNTracks) delete fVolNTracks;
269 if(fhEmpty) delete fhEmpty;
270 if(fhistVolNptsUsed) delete fhistVolNptsUsed;
271 if(fhistVolUsed) delete fhistVolUsed;
272 if(fSigmaVolZ) delete fSigmaVolZ;
273 if(fpTrackVolIDs) delete fpTrackVolIDs;
274 if(fVolVolids) delete fVolVolids;
275 if(fVolUsed) delete fVolUsed;
276 if(fClonesArray) delete fClonesArray;
278 SetFileNameTrackPoints("");
279 SetFileNameGeometry("");
283 //____________________________________________________________________________
284 void AliITSResidualsAnalysis::InitHistograms(const TArrayI *volIDs)
287 // Method that sets and creates the required hisstrograms
288 // with the correct binning (it dos not fill them)
292 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
294 TString histnameRPHI="HistRPHI_volID_",aux;
295 TString histnameZ="HistZ_volID_";
296 TString histnameX="HistX_volID_";
297 TString histnameGlob="HistGlob_volID_";
298 TString histnameCorrVol="HistCorrVol_volID";
299 TString histnamePullRPHI="HistPullRPHI_volID_";
300 TString histnamePullZ="HistPullZ_volID_";
302 TString histnameDirPhi="HistTrackDirPhi_volID_";
303 TString histnameDirLambda="HistTrackDirLambda_volID_";
304 TString histnameDirLambda2="HistTrackDirLambda2_volID_";
305 TString histnameDirAlpha="HistTrackDirAlpha_volID_";
306 TString histnameDir="HistTrackDir_volID_";
308 fnHist=volIDs->GetSize();
309 fVolResHistRPHI=new TH1F*[fnHist];
310 fResHistGlob=new TH1F*[fnHist];
311 fResHistZ=new TH1F*[fnHist];
312 fResHistX=new TH1F*[fnHist];
313 fResHistXLocsddL=new TH1F*[fnHist];
314 fResHistXLocsddR=new TH1F*[fnHist];
315 fHistCoordGlobY=new TH1F*[fnHist];
316 fPullHistRPHI=new TH1F*[fnHist];
317 fPullHistZ=new TH1F*[fnHist];
318 fhistCorrVol=new TH2F*[fnHist];
320 fTrackDirPhi=new TH1F*[fnHist];
321 fTrackDirLambda=new TH1F*[fnHist];
322 fTrackDirLambda2=new TH1F*[fnHist];
323 fTrackDirAlpha=new TH1F*[fnHist];
325 fTrackDirPhiAll=new TH1F("fTrackDirPhiAll","fTrackDirPhiAll",100,-180,180);
326 fTrackDirLambdaAll=new TH1F("fTrackDirLambdaAll","fTrackDirLambdaAll",100,-180,180);
327 fTrackDirLambda2All=new TH1F("fTrackDirLambda2All","fTrackDirLambda2All",100,0,180);
328 fTrackDirAlphaAll=new TH1F("fTrackDirAlphaAll","fTrackDirAlphaAll",100,-180,180);
330 fTrackDirAll=new TH2F("fTrackDirAll","Hist with trakcs directions",100,-180,180,100,-180,180);
331 fTrackDir2All=new TH2F("fTrackDir2All","Hist with trakcs directions",100,-180,180,100,-180,180);
332 fTrackDirXZAll=new TH2F("fTrackDirXZAll","Hist with trakcs directions from XZ ",100,-3.,3.,100,-3.,3.);
334 fTrackDir=new TH2F*[fnHist];
337 Float_t **binningZPhi;
341 binningZPhi=CheckSingleLayer(volIDs);
342 fvolidsToBin=new Int_t*[fnPhi*fnZ];
343 binningphi=binningZPhi[0];
344 binningz=binningZPhi[1];
345 binning=SetBinning(volIDs,binningphi,binningz);
347 if(binning){ //ONLY FOR A SINGLE LAYER!
348 fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
349 fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",fnPhi,binningphi,fnZ,binningz);
350 fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",fnPhi,binningphi,fnZ,binningz);
351 fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",fnPhi,binningphi,fnZ,binningz);
352 fhEmpty=new TH2F("fhEmpty","Hist for getting binning",fnPhi,binningphi,fnZ,binningz);
353 fVolNTracks->SetXTitle("Volume #phi");
354 fVolNTracks->SetYTitle("Volume z ");
355 fhistVolNptsUsed->SetXTitle("Volume #phi");
356 fhistVolNptsUsed->SetYTitle("Volume z ");
357 fhistVolUsed->SetXTitle("Volume #phi");
358 fhistVolUsed->SetYTitle("Volume z ");
359 fSigmaVolZ->SetXTitle("Volume #phi");
360 fSigmaVolZ->SetYTitle("Volume z ");
362 fVolNTracks=new TH2F("fVolNTracks","Hist with N tracks passing through a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
363 fhistVolNptsUsed=new TH2F("fhistVolNptsUsed","Hist with N points used for given module==(r,phi) ",50,-3.2,3.2,100,-80,80);
364 fhistVolUsed=new TH2F("fhistVolUsed","Hist with N modules used for a given module==(r,phi) zone",50,-3.2,3.2,100,-80,80);
365 fSigmaVolZ=new TH2F("fSigmaVolZ","Hist with Sigma of Residual distribution for each module",50,-3.2,3.2,100,-80,80);
366 fhEmpty=new TH2F("fhEmpty","Hist for getting binning",50,-3.2,3.2,100,-80,80);
367 fVolNTracks->SetXTitle("Volume #phi");
368 fVolNTracks->SetYTitle("Volume z ");
369 fhistVolNptsUsed->SetXTitle("Volume #phi");
370 fhistVolNptsUsed->SetYTitle("Volume z ");
371 fhistVolUsed->SetXTitle("Volume #phi");
372 fhistVolUsed->SetYTitle("Volume z ");
373 fSigmaVolZ->SetXTitle("Volume #phi");
374 fSigmaVolZ->SetYTitle("Volume z ");
377 fpTrackVolIDs=new TArrayI(fnHist);
378 fVolUsed=new TArrayI*[fnHist];
379 fVolVolids=new TArrayI*[fnHist];
380 fLastVolVolid=new Int_t[fnHist];
382 for (Int_t nhist=0;nhist<fnHist;nhist++){
383 fpTrackVolIDs->AddAt(volIDs->At(nhist),nhist);
385 aux+=volIDs->At(nhist);
386 fVolResHistRPHI[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
387 fVolResHistRPHI[nhist]->SetName(aux.Data());
388 fVolResHistRPHI[nhist]->SetTitle(aux.Data());
391 aux+=volIDs->At(nhist);
392 fResHistZ[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
393 fResHistZ[nhist]->SetName(aux.Data());
394 fResHistZ[nhist]->SetTitle(aux.Data());
397 aux+=volIDs->At(nhist);
398 fResHistX[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
399 fResHistX[nhist]->SetName(aux.Data());
400 fResHistX[nhist]->SetTitle(aux.Data());
403 aux+=volIDs->At(nhist);
404 aux.Append("LocalSDDLeft");
405 fResHistXLocsddL[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
406 fResHistXLocsddL[nhist]->SetName(aux.Data());
407 fResHistXLocsddL[nhist]->SetTitle(aux.Data());
410 aux+=volIDs->At(nhist);
411 aux.Append("LocalSDDRight");
412 fResHistXLocsddR[nhist]=new TH1F("histname","histname",20000,-5.0,5.0);
413 fResHistXLocsddR[nhist]->SetName(aux.Data());
414 fResHistXLocsddR[nhist]->SetTitle(aux.Data());
416 aux="fHistCoordGlobY";
417 fHistCoordGlobY[nhist]=new TH1F("histname","histname",24000,-30.,30.);
418 fHistCoordGlobY[nhist]->SetName(aux.Data());
419 fHistCoordGlobY[nhist]->SetTitle(aux.Data());
421 aux=histnamePullRPHI;
422 aux+=volIDs->At(nhist);
423 fPullHistRPHI[nhist]=new TH1F("histname","histname",100,-7.,7.);
424 fPullHistRPHI[nhist]->SetName(aux.Data());
425 fPullHistRPHI[nhist]->SetTitle(aux.Data());
428 aux+=volIDs->At(nhist);
429 fPullHistZ[nhist]=new TH1F("histname","histname",100,-7.,7.);
430 fPullHistZ[nhist]->SetName(aux.Data());
431 fPullHistZ[nhist]->SetTitle(aux.Data());
434 aux+=volIDs->At(nhist);
435 fTrackDirPhi[nhist]=new TH1F("histname","histname",100,-180,180);
436 fTrackDirPhi[nhist]->SetName(aux.Data());
437 fTrackDirPhi[nhist]->SetTitle(aux.Data());
439 aux=histnameDirLambda;
440 aux+=volIDs->At(nhist);
441 fTrackDirLambda[nhist]=new TH1F("histname","histname",100,0,180);
442 fTrackDirLambda[nhist]->SetName(aux.Data());
443 fTrackDirLambda[nhist]->SetTitle(aux.Data());
445 aux=histnameDirLambda2;
446 aux+=volIDs->At(nhist);
447 fTrackDirLambda2[nhist]=new TH1F("histname","histname",100,0,180);
448 fTrackDirLambda2[nhist]->SetName(aux.Data());
449 fTrackDirLambda2[nhist]->SetTitle(aux.Data());
451 aux=histnameDirAlpha;
452 aux+=volIDs->At(nhist);
453 fTrackDirAlpha[nhist]=new TH1F("histname","histname",100,-180,180);
454 fTrackDirAlpha[nhist]->SetName(aux.Data());
455 fTrackDirAlpha[nhist]->SetTitle(aux.Data());
458 aux+=volIDs->At(nhist);
459 fTrackDir[nhist]=new TH2F("histname","histname",100,-90.,90.,100,-180.,180.);
460 fTrackDir[nhist]->SetName(aux.Data());
461 fTrackDir[nhist]->SetTitle(aux.Data());
464 aux+=volIDs->At(nhist);
465 fResHistGlob[nhist]=new TH1F("histname","histname",400,-0.08,0.08);
466 fResHistGlob[nhist]->SetName(aux.Data());
467 fResHistGlob[nhist]->SetTitle(aux.Data());
470 aux+=volIDs->At(nhist);
471 fhistCorrVol[nhist]=new TH2F("histname","histname",50,-3.2,3.2,100,-80,80);
474 fhistCorrVol[nhist]->SetName(aux.Data());
475 fhistCorrVol[nhist]->SetTitle(aux.Data());
476 fhistCorrVol[nhist]->SetXTitle("Volume #varphi");
477 fhistCorrVol[nhist]->SetYTitle("Volume z ");
478 fVolVolids[nhist]=new TArrayI(100);
479 fVolUsed[nhist]=new TArrayI(1000);
480 fLastVolVolid[nhist]=0;
488 //____________________________________________________________________________
489 void AliITSResidualsAnalysis::ListVolUsed(TTree *pointsTree,TArrayI ***arrayIndex,Int_t **lastIndex)
492 // This is copied from AliAlignmentClass::LoadPoints() method
494 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
495 Int_t volIDalignable,volIDpoint,iModule;
497 AliTrackPointArray* array = 0;
498 pointsTree->SetBranchAddress("SP", &array);
501 for(Int_t ivol=0;ivol<fnHist;ivol++){
503 volIDalignable=fpTrackVolIDs->At(ivol);
504 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volIDalignable,iModule);
506 Int_t nArraysId = lastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
507 TArrayI *index = arrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
508 for (Int_t iArrayId = 0;iArrayId < nArraysId; iArrayId++) {
511 Int_t entry = (*index)[iArrayId];
513 pointsTree->GetEvent(entry);
515 AliWarning("Wrong space point array index!");
519 // Get the space-point array
520 Int_t modnum,nPoints = array->GetNPoints();
522 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
525 array->GetPoint(p,iPoint);
527 AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
528 // check if the layer id is valid
529 if ((layer < AliGeomManager::kFirstLayer) ||
530 (layer >= AliGeomManager::kLastLayer)) {
531 AliError(Form("Layer index is invalid: %d (%d -> %d) !",
532 layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
537 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
539 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
540 layer,modnum,AliGeomManager::LayerSize(layer)));
543 if (layer > AliGeomManager::kSSD2) continue; // ITS only
546 volIDpoint=(Int_t)p.GetVolumeID();
547 if (volIDpoint==volIDalignable) continue;
548 Int_t size = fVolVolids[ivol]->GetSize();
549 // If needed allocate new size
550 if (fLastVolVolid[ivol]>=size){// Warning: fLAST[NHIST] is useless
551 fVolVolids[ivol]->Set(size + 1000);
555 fVolVolids[ivol]->AddAt(volIDpoint,fLastVolVolid[ivol]);
556 fLastVolVolid[ivol]++;
557 Bool_t usedVol=kFALSE;
558 for(Int_t used=0;used<lastused;used++){
559 if(fVolUsed[ivol]->At(used)==volIDpoint){
567 size = fVolUsed[ivol]->GetSize();
568 // If needed allocate new size
569 if (lastused>= size){
570 fVolUsed[ivol]->Set(size + 1000);
572 fVolUsed[ivol]->AddAt(volIDpoint,lastused);
576 FillVolumeCorrelationHists(ivol,volIDalignable,volIDpoint,usedVol);
585 //____________________________________________________________________________
586 void AliITSResidualsAnalysis::FillVolumeCorrelationHists(Int_t ivol,Int_t volIDalignable,Int_t volIDpoint,Bool_t usedVol) const
589 // Fill the histograms with the correlations between volumes
593 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
594 Double_t translGlobal[3];
597 // const char *symname,*volpath;
599 TGeoPhysicalNode *pn;
600 TGeoHMatrix *globMatrix;
603 symname = AliGeomManager::SymName(volIDalignable);
604 pne = gGeoManager->GetAlignableEntry(symname);
605 volpath=pne->GetTitle();
606 pn=gGeoManager->MakePhysicalNode(volpath);
607 globMatrix=pn->GetMatrix();
610 AliGeomManager::GetOrigTranslation(volIDalignable,translGlobal);
611 // radius=TMath::Sqrt(transGlobal[0]*transGlobal[0]+transGlobal[1]*transGlobal[1]);
612 phi=TMath::ATan2(translGlobal[1],translGlobal[0]);
613 fhistVolNptsUsed->Fill(phi,translGlobal[2]);
615 fhistVolUsed->Fill(phi,translGlobal[2]);
618 /* symname = AliGeomManager::SymName(volIDpoint);
619 pne = gGeoManager->GetAlignableEntry(symname);
620 volpath=pne->GetTitle();
621 pn=gGeoManager->MakePhysicalNode(volpath);
622 globMatrix=pn->GetMatrix();
623 transGlobal=globMatrix->GetTranslation();
625 AliGeomManager::GetOrigTranslation(volIDpoint,translGlobal);
626 // radius=TMath::Sqrt(transGlobal[0]*transGlobal[0]+transGlobal[1]*transGlobal[1]);
627 phi=TMath::ATan2(translGlobal[1],translGlobal[0]);
629 fhistCorrVol[ivol]->Fill(phi,translGlobal[2]);
634 //____________________________________________________________________________
635 void AliITSResidualsAnalysis::FillResidualsH(AliTrackPointArray *points,
636 AliTrackPointArray *pTrack) const
639 // Method that fills the histograms with the residuals
643 Float_t xyz[3],xyz2[3];
644 Double_t xyzD[3],xyz2D[3];
645 Double_t loc[3],loc2[3];
647 Float_t resRPHI,resGlob,resZ,resX;
649 Double_t pullrphi,sign,phi;
652 for(Int_t ipoint=0;ipoint<points->GetNPoints();ipoint++){
654 //pTrack->GetPoint(pTr,ipoint);
655 points->GetPoint(p,ipoint);
656 volIDpoint=(Int_t)p.GetVolumeID();
659 pTrack->GetPoint(pTr,ipoint);
662 for(Int_t i=0;i<3;i++){
667 phi = TMath::ATan2(xyz[1],xyz[0]);//<-watch out: phi of the pPoints!
672 resRPHI=TMath::Sqrt((xyz2[0]-xyz[0])*(xyz2[0]-xyz[0])+(xyz2[1]-xyz[1])*(xyz2[1]-xyz[1]));
674 sign=TMath::ATan2(xyz2[1],xyz2[0])-TMath::ATan2(xyz[1],xyz[0]);
676 sign=sign/TMath::Abs(sign);
677 resRPHI=resRPHI*sign;
685 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]));
687 for(Int_t ivolIDs=0;ivolIDs<fpTrackVolIDs->GetSize();ivolIDs++){
688 if(volIDpoint==fpTrackVolIDs->At(ivolIDs)){
690 fVolResHistRPHI[ivolIDs]->Fill(resRPHI);
691 fResHistZ[ivolIDs]->Fill(resZ);
692 fResHistX[ivolIDs]->Fill(resX);
693 fHistCoordGlobY[ivolIDs]->Fill(xyz[1]);
695 Int_t modIndex = -1; // SDD Section
696 if(AliGeomManager::VolUIDToLayer(volIDpoint)==3) modIndex=volIDpoint-6144+240;
697 if(AliGeomManager::VolUIDToLayer(volIDpoint)==4) modIndex=volIDpoint-8192+240+84;
699 AliITSgeomTGeo::GlobalToLocal(modIndex,xyzD,loc); // error here!?
700 AliITSgeomTGeo::GlobalToLocal(modIndex,xyz2D,loc2);
701 Float_t rexloc=loc2[0]-loc[0];
702 //cout<<"Residual: "<<volIDpoint<<" "<<loc[0]<<" -> "<<rexloc<<endl;
704 fResHistXLocsddR[ivolIDs]->Fill(rexloc);
706 fResHistXLocsddL[ivolIDs]->Fill(rexloc);
709 fResHistGlob[ivolIDs]->Fill(resGlob);
711 fTrackDirPhiAll->Fill(phi);
712 fTrackDirPhi[ivolIDs]->Fill(phi);
716 Float_t globalPhi,globalZ;
717 if(kTRUE||(fvolidsToBin[ivolIDs][0]!=volIDpoint)){
718 binphi=GetBinPhiZ((Int_t)volIDpoint,&binz);
721 // This in the case of alignment of one entire layer
722 // (fnHIst=layersize) may reduce iterations:
723 // remind of that fsingleLayer->fnHista<layerSize
724 binphi=fvolidsToBin[ivolIDs][1];
725 binz=fvolidsToBin[ivolIDs][2];
727 globalPhi=fCoordToBinTable[binphi][binz][0];
728 globalZ=fCoordToBinTable[binphi][binz][1];
730 fVolNTracks->Fill(globalPhi,globalZ);
732 else fVolNTracks->Fill(TMath::ATan2(xyz[1],xyz[0]),xyz[2]);
738 //____________________________________________________________________________
739 Bool_t AliITSResidualsAnalysis::SaveHists(Int_t minNpoints, TString outname) const
742 // Saves the histograms into a tree and saves the tree into a file
747 TFile *hFile=new TFile(outname.Data(),"RECREATE","File containing the Residuals Tree");
749 // TTree with the residuals
750 TTree *analysisTree=new TTree("analysisTree","Tree with the residuals");
752 // Declares Variables to be stored into the TTree
753 TF1 *gauss=new TF1("gauss","gaus",-10.,10.);
754 Int_t volID,entries,nHistAnalyzed=0;
755 Double_t meanResRPHI,meanResZ,meanResX,rmsResRPHI,rmsResZ,rmsResX,coordVol[3],x,y,z;
756 TH1F *histRPHI = new TH1F();
757 TH1F *histZ = new TH1F();
758 TH1F *histX = new TH1F();
759 TH1F *histXLocsddL = new TH1F();
760 TH1F *histXLocsddR = new TH1F();
761 TH1F *histCoordGlobY = new TH1F();
762 // Note: 0 = RPHI, 1 = Z
765 // Branching the TTree
766 analysisTree->Branch("volID",&volID,"volID/I");
767 analysisTree->Branch("x",&x,"x/D");
768 analysisTree->Branch("y",&y,"y/D");
769 analysisTree->Branch("z",&z,"z/D");
770 analysisTree->Branch("meanResRPHI",&meanResRPHI,"meanResRPHI/D");
771 analysisTree->Branch("meanResZ",&meanResZ,"meanResZ/D");
772 analysisTree->Branch("meanResX",&meanResX,"meanResX/D");
773 analysisTree->Branch("rmsResRPHI",&rmsResRPHI,"rmsResRPHI/D");
774 analysisTree->Branch("rmsResZ",&rmsResZ,"rmsResZ/D");
776 analysisTree->Branch("histRPHI","TH1F",&histRPHI,128000,0);
777 analysisTree->Branch("histZ","TH1F",&histZ,128000,0);
778 analysisTree->Branch("histX","TH1F",&histX,128000,0);
779 analysisTree->Branch("histXLocsddL","TH1F",&histXLocsddL,128000,0);
780 analysisTree->Branch("histXLocsddR","TH1F",&histXLocsddR,128000,0);
781 analysisTree->Branch("histCoordGlobY","TH1F",&histCoordGlobY,128000,0);
785 for(Int_t j=0;j<fnHist;j++){
787 volID=fpTrackVolIDs->At(j);
788 AliGeomManager::GetTranslation(volID,coordVol);
793 entries=(Int_t)(fResHistGlob[j]->GetEntries());
796 if(entries>=minNpoints){
800 //entries=(Int_t)fVolResHistRPHI[j]->GetEntries();
803 histRPHI=fVolResHistRPHI[j];
804 rmsResRPHI=fVolResHistRPHI[j]->GetRMS();
807 gauss->SetRange(-3*rmsResRPHI,3*rmsResRPHI);
808 fVolResHistRPHI[j]->Fit("gauss","QRN");
809 meanResRPHI=gauss->GetParameter(1);
811 meanResRPHI=fVolResHistRPHI[j]->GetMean();
816 rmsResZ=fResHistZ[j]->GetRMS();
819 gauss->SetRange(-3*rmsResZ,3*rmsResZ);
820 fResHistZ[j]->Fit("gauss","QRN");
821 meanResZ=gauss->GetParameter(1);
823 meanResZ=fResHistZ[j]->GetMean();
828 rmsResX=fResHistX[j]->GetRMS();
831 gauss->SetRange(-3*rmsResX,3*rmsResX);
832 fResHistX[j]->Fit("gauss","QRN");
833 meanResX=gauss->GetParameter(1);
835 meanResX=fResHistX[j]->GetMean();
838 histXLocsddL=fResHistXLocsddL[j];
839 histXLocsddR=fResHistXLocsddR[j];
840 histCoordGlobY=fHistCoordGlobY[j];
842 analysisTree->Fill();
846 //entries=(Int_t)fVolResHistRPHI[j]->GetEntries();
849 histRPHI=fVolResHistRPHI[j];
862 histXLocsddL=fResHistXLocsddL[j];
863 histXLocsddR=fResHistXLocsddR[j];
864 histCoordGlobY=fHistCoordGlobY[j];
866 analysisTree->Fill();
873 cout<<"-> Modules Analyzed: "<<nHistAnalyzed<<endl;
874 cout<<" With "<<blimps<<" events"<<endl;
878 analysisTree->Write();
879 fVolNTracks->Write();
882 //TCanvas *color = new TCanvas("color","fhistVolUsed",800,600);
883 //fhistVolUsed->DrawCopy("COLZ");
885 fhistVolUsed->Write();
886 /* fTrackDirPhiAll->Write();
887 fTrackDirLambdaAll->Write();
888 fTrackDirLambda2All->Write();
889 fTrackDirAlphaAll->Write();
890 fTrackDirAll->Write();
891 fTrackDir2All->Write();
892 fTrackDirXZAll->Write();
894 fhistVolNptsUsed->Write();
895 hFile->mkdir("CorrVol");
896 hFile->cd("CorrVol");
897 for(Int_t corr=0;corr<fnHist;corr++)fhistCorrVol[corr]->Write();
900 // fhistVolNptsUsed->Write();
909 //____________________________________________________________________________
910 void AliITSResidualsAnalysis::DrawHists() const
913 // Draws the histograms of the residuals and of the number of tracks
917 for(Int_t canv=0;canv<fnHist;canv++){
918 cname="canv Residuals";
920 TCanvas *c=new TCanvas(cname.Data(),cname.Data(),700,700);
923 fVolResHistRPHI[canv]->Draw();
925 fResHistZ[canv]->Draw();
927 fResHistGlob[canv]->Draw();
929 cname="canv NVolTracks";
931 TCanvas *c2=new TCanvas(cname.Data(),cname.Data(),700,700);
938 //____________________________________________________________________________
939 Float_t** AliITSResidualsAnalysis::CheckSingleLayer(const TArrayI *volids)
942 // Checks if volumes array is a single (ITS) layer or not
945 Float_t **binningzphi=new Float_t*[2];
947 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volids->At(0),iModule);
948 //Check that one single Layer is going to be aligned
949 for(Int_t nvol=0;nvol<volids->GetSize();nvol++){
950 if(iLayer != AliGeomManager::VolUIDToLayer((UShort_t)volids->At(nvol),iModule)){
951 printf("Wrong Layer! \n %d , %d , %d ,%d \n",(Int_t)AliGeomManager::VolUIDToLayer((UShort_t)volids->At(nvol),iModule),nvol,volids->GetSize(),iModule);
957 //Bool_t used=kFALSE;
959 case AliGeomManager::kSPD1:{
960 fnPhi=kPhiSPD1;//kPhiSPD1;
962 binningzphi[0]=new Float_t[kPhiSPD1+1];
963 binningzphi[1]=new Float_t[kZSPD1+1];
964 fCoordToBinTable=new Double_t**[kPhiSPD1];
965 for(Int_t j=0;j<fnPhi;j++){
966 fCoordToBinTable[j]=new Double_t*[kZSPD1];
969 case AliGeomManager::kSPD2:{
970 fnPhi=kPhiSPD2;//kPhiSPD1;
972 binningzphi[0]=new Float_t[kPhiSPD2+1];
973 binningzphi[1]=new Float_t[kZSPD2+1];
974 fCoordToBinTable=new Double_t**[kPhiSPD2];
975 for(Int_t j=0;j<fnPhi;j++){
976 fCoordToBinTable[j]=new Double_t*[kZSPD2];
979 }; break; case AliGeomManager::kSDD1:{
980 fnPhi=kPhiSDD1;//kPhiSPD1;
982 binningzphi[0]=new Float_t[kPhiSDD1+1];
983 binningzphi[1]=new Float_t[kZSDD1+1];
984 fCoordToBinTable=new Double_t**[kPhiSDD1];
985 for(Int_t j=0;j<fnPhi;j++){
986 fCoordToBinTable[j]=new Double_t*[kZSDD1];
988 }; break; case AliGeomManager::kSDD2:{
989 fnPhi=kPhiSDD2;//kPhiSPD1;
991 binningzphi[0]=new Float_t[kPhiSDD2+1];
992 binningzphi[1]=new Float_t[kZSDD2+1];
993 fCoordToBinTable=new Double_t**[kPhiSDD2];
994 for(Int_t j=0;j<fnPhi;j++){
995 fCoordToBinTable[j]=new Double_t*[kZSDD2];
997 }; break; case AliGeomManager::kSSD1:{
998 fnPhi=kPhiSSD1;//kPhiSPD1;
1000 binningzphi[0]=new Float_t[kPhiSSD1+1];
1001 binningzphi[1]=new Float_t[kZSSD1+1];
1002 fCoordToBinTable=new Double_t**[kPhiSSD1];
1003 for(Int_t j=0;j<fnPhi;j++){
1004 fCoordToBinTable[j]=new Double_t*[kZSSD1];
1006 }; break; case AliGeomManager::kSSD2:{
1007 fnPhi=kPhiSSD2;//kPhiSPD1;
1008 fnZ=kZSSD2;//nZSPD1;
1009 binningzphi[0]=new Float_t[kPhiSSD2+1];
1010 binningzphi[1]=new Float_t[kZSSD2+1];
1011 fCoordToBinTable=new Double_t**[kPhiSSD2];
1012 for(Int_t j=0;j<fnPhi;j++){
1013 fCoordToBinTable[j]=new Double_t*[kZSSD2];
1018 printf("Wrong Layer Label! \n");
1019 fsingleLayer=kFALSE;
1020 //////////////////////
1023 binningzphi[0]=new Float_t[1];
1024 binningzphi[1]=new Float_t[1];
1025 fCoordToBinTable=new Double_t**[1];
1026 for(Int_t j=0;j<fnPhi;j++){
1027 fCoordToBinTable[j]=new Double_t*[1];
1030 /////////////////////
1038 //____________________________________________________________________________
1039 Bool_t AliITSResidualsAnalysis::SetBinning(const TArrayI *volids,Float_t *phiBin,Float_t *zBin)
1042 // Sets the correct binning
1045 if(!fsingleLayer)return kFALSE;
1046 // const char *volpath,*symname;
1048 Int_t *orderArrayPhi,*orderArrayZ;
1050 Double_t *phiArray,*zArray,*phiArrayOrdered,*zArrayOrdered;
1051 Double_t translGlobal[3];
1052 Double_t lastPhimin=-10;
1053 Double_t lastZmin=-99;
1055 /* TGeoPNEntry *pne;
1056 TGeoPhysicalNode *pn;
1057 TGeoHMatrix *globMatrix;*/
1061 orderPhiZ=new Int_t**[fnPhi];
1062 phiArray=new Double_t[fnPhi];//phiBin[nModulesPhi+1];
1063 zArray=new Double_t[fnZ];//zBin[nModulesZ+1];
1064 phiArrayOrdered=new Double_t[fnPhi];
1065 zArrayOrdered=new Double_t[fnZ];
1066 orderArrayPhi=new Int_t[fnPhi];
1067 orderArrayZ=new Int_t[fnZ];
1068 for(Int_t k=0;k<fnZ;k++){
1072 for(Int_t k=0;k<fnPhi;k++){
1075 orderPhiZ[k]=new Int_t*[fnZ];
1079 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer((UShort_t)volids->At(0),iModule);
1081 Int_t lastPhi=0,lastZ=0;
1082 for(iModule=0;iModule<AliGeomManager::LayerSize(iLayer);iModule++){
1083 fvolidsToBin[iModule]=new Int_t[3];
1084 volID=AliGeomManager::LayerToVolUID(iLayer,iModule);
1085 fvolidsToBin[iModule][0]=volID;
1086 /* symname = AliGeomManager::SymName(volID);
1087 pne = gGeoManager->GetAlignableEntry(symname);
1088 volpath=pne->GetTitle();
1089 pn=gGeoManager->MakePhysicalNode(volpath);
1090 globMatrix=pn->GetMatrix();
1091 translGlobal=globMatrix->GetTranslation();
1094 AliGeomManager::GetOrigTranslation(volID,translGlobal);
1096 for(Int_t j=0;j<lastPhi;j++){
1098 if(TMath::Abs(phiArray[j]-TMath::ATan2(translGlobal[1],translGlobal[0]))<2*TMath::Pi()/(10*fnPhi)){//10 is a safety factor but....
1099 fvolidsToBin[iModule][1]=j;
1105 phiArray[lastPhi]=TMath::ATan2(translGlobal[1],translGlobal[0]);
1106 fvolidsToBin[iModule][1]=lastPhi;
1107 if(phiArray[lastPhi]<lastPhimin)lastPhimin=phiArray[lastPhi];
1110 printf("Wrong Phi! \n");
1114 for(Int_t j=0;j<lastZ;j++){
1116 if(TMath::Abs(zArray[j]-translGlobal[2])<0.1){
1117 fvolidsToBin[iModule][2]=j;
1123 fvolidsToBin[iModule][2]=lastZ;
1124 zArray[lastZ]=translGlobal[2];
1125 if(zArray[lastZ]<lastZmin)lastZmin=zArray[lastZ];
1128 printf("Wrong Z! \n");
1134 //ORDERING THE ARRAY OF PHI AND Z VALUES
1135 for(Int_t order=0;order<fnPhi;order++){
1136 for(Int_t j=0;j<fnPhi;j++){
1137 if((j!=order)&&(phiArray[order]>phiArray[j]))orderArrayPhi[order]++;
1141 for(Int_t order=0;order<fnPhi;order++){
1142 for(Int_t j=0;j<fnPhi;j++){
1143 if(orderArrayPhi[j]==order){
1144 phiArrayOrdered[order]=phiArray[j];
1151 for(Int_t order=0;order<fnZ;order++){
1152 for(Int_t j=0;j<fnZ;j++){
1153 if((j!=order)&&(zArray[order]>zArray[j]))orderArrayZ[order]++;
1158 for(Int_t order=0;order<fnZ;order++){
1159 for(Int_t j=0;j<fnZ;j++){
1160 if(orderArrayZ[j]==order){
1161 zArrayOrdered[order]=zArray[j];
1168 //Filling the fCoordToBinTable
1169 for(Int_t j=0;j<fnPhi;j++){
1170 for(Int_t i=0;i<fnZ;i++){
1171 orderPhiZ[j][i]=new Int_t[2];
1172 orderPhiZ[j][i][0]=orderArrayPhi[j];
1173 orderPhiZ[j][i][1]=orderArrayZ[i];
1178 for(Int_t j=0;j<fnPhi;j++){
1179 for(Int_t i=0;i<fnZ;i++){
1180 fCoordToBinTable[j][i]=new Double_t[2];
1181 fCoordToBinTable[j][i][0]=phiArrayOrdered[j];
1182 fCoordToBinTable[j][i][1]=zArrayOrdered[i];
1183 printf("Now (binphi,binz)= %d, %d e (phi,z)=%f,%f \n",j,i,fCoordToBinTable[j][i][0],fCoordToBinTable[j][i][1]);
1187 for(iModule=0;iModule<fnPhi*fnZ;iModule++){
1188 istar=fvolidsToBin[iModule][1];
1189 jstar=fvolidsToBin[iModule][2];
1190 fvolidsToBin[iModule][1]=orderPhiZ[istar][jstar][0];
1191 fvolidsToBin[iModule][2]=orderPhiZ[istar][jstar][1];
1195 //now constructing the binning
1196 for(Int_t iModPhi=0;iModPhi<fnPhi-1;iModPhi++){
1197 phiBin[iModPhi+1]=(Float_t)phiArrayOrdered[iModPhi]+0.5*(phiArrayOrdered[iModPhi+1]-phiArrayOrdered[iModPhi]);
1200 phiBin[0]=(Float_t)phiArrayOrdered[0]-(phiArrayOrdered[1]-phiArrayOrdered[0])/2.;
1202 phiBin[fnPhi]=(Float_t)phiArrayOrdered[fnPhi-1]+(phiArrayOrdered[fnPhi-1]-phiArrayOrdered[fnPhi-2])/2.;
1203 for(Int_t iModPhi=0;iModPhi<fnPhi+1;iModPhi++){
1204 printf("Mean Phi mod %d su %d: %f \n",iModPhi+1,fnPhi,phiBin[iModPhi]);
1207 for(Int_t iModZ=0;iModZ<fnZ-1;iModZ++){
1208 zBin[iModZ+1]=(Float_t)zArrayOrdered[iModZ]+0.5*(zArrayOrdered[iModZ+1]-zArrayOrdered[iModZ]);
1210 zBin[0]=(Float_t)zArrayOrdered[0]-0.5*(zArrayOrdered[1]-zArrayOrdered[0]);
1211 zBin[fnZ]=(Float_t)zArrayOrdered[fnZ-1]+0.5*(zArrayOrdered[1]-zArrayOrdered[0]);
1214 for(Int_t iModPhi=0;iModPhi<fnZ+1;iModPhi++){
1215 printf("Mean Z mod %d su %d: %f \n",iModPhi+1,fnZ,zBin[iModPhi]);
1221 //____________________________________________________________________________
1222 Int_t AliITSResidualsAnalysis::GetBinPhiZ(const Int_t volID,Int_t *binz) const
1225 // Returns the correct Phi-Z bin
1229 printf("No Single Layer reAlignment! \n");
1233 for(Int_t j=0;j<fnPhi*fnZ;j++){
1235 printf("Wrong volume UID introduced! fnHist: %d volID: %d iter: %d \n",fnHist,volID,j);
1238 if(fvolidsToBin[j][0]==volID){
1240 *binz=fvolidsToBin[j][2];
1241 return fvolidsToBin[j][1];
1249 //___________________________________________________________________________
1250 Int_t AliITSResidualsAnalysis::WhichSector(Int_t module) const
1253 // This method returns the number of the SPD Sector
1254 // to which belongs the module (Sectors 0-9)
1256 //--->cSect = 0 <---
1280 || module==4111) return 0;
1282 //--->cSect = 1 <---
1306 || module==4127) return 1;
1308 //--->cSect = 2 <---
1332 || module==4143) return 2;
1334 //--->cSect = 3 <---
1358 || module==4159) return 3;
1360 //--->cSect = 4 <---
1384 || module==4175) return 4;
1386 //--->cSect = 5 <---
1410 || module==4191) return 5;
1412 //--->cSect = 6 <---
1436 || module==4207) return 6;
1438 //--->cSect = 7 <---
1462 || module==4223) return 7;
1464 //--->cSect = 8 <---
1488 || module==4239) return 8;
1490 //--->cSect = 9 <---
1514 || module==4255) return 9;
1516 //printf("Module not belonging to SPD, sorry!");
1521 //____________________________________________________________________________
1522 TArrayI* AliITSResidualsAnalysis::GetSPDSectorsVolids(Int_t sectors[10]) const
1525 // This method gets the volID Array for the chosen sectors.
1526 // You have to pass an array with a 1 for each selected sector.
1527 // i.e. sectors[10] = {1,1,0,0,0,0,0,0,1,0} -> Sector 0, 1, 9 selected.
1533 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1535 for(Int_t co=0;co<10;co++){ //counts the number of sectors chosen
1536 if(sectors[co]==1) nSect++;
1539 if(nSect<1){ //if no sector chosen -> exit
1540 Printf("Error! No Sector/s Selected!");
1544 TArrayI *volIDs = new TArrayI(nSect*24);
1546 if(sectors[0]==1){ //--->cSect = 0 <---
1547 volIDs->AddAt(2048,iModule); iModule++;
1548 volIDs->AddAt(2049,iModule); iModule++;
1549 volIDs->AddAt(2050,iModule); iModule++;
1550 volIDs->AddAt(2051,iModule); iModule++;
1551 volIDs->AddAt(2052,iModule); iModule++;
1552 volIDs->AddAt(2053,iModule); iModule++;
1553 volIDs->AddAt(2054,iModule); iModule++;
1554 volIDs->AddAt(2055,iModule); iModule++;
1555 volIDs->AddAt(4096,iModule); iModule++;
1556 volIDs->AddAt(4097,iModule); iModule++;
1557 volIDs->AddAt(4098,iModule); iModule++;
1558 volIDs->AddAt(4099,iModule); iModule++;
1559 volIDs->AddAt(4100,iModule); iModule++;
1560 volIDs->AddAt(4101,iModule); iModule++;
1561 volIDs->AddAt(4102,iModule); iModule++;
1562 volIDs->AddAt(4103,iModule); iModule++;
1563 volIDs->AddAt(4104,iModule); iModule++;
1564 volIDs->AddAt(4105,iModule); iModule++;
1565 volIDs->AddAt(4106,iModule); iModule++;
1566 volIDs->AddAt(4107,iModule); iModule++;
1567 volIDs->AddAt(4108,iModule); iModule++;
1568 volIDs->AddAt(4109,iModule); iModule++;
1569 volIDs->AddAt(4110,iModule); iModule++;
1570 volIDs->AddAt(4111,iModule); iModule++;
1572 if(sectors[1]==1){ //--->cSect = 1 <//---
1573 volIDs->AddAt(2056,iModule); iModule++;
1574 volIDs->AddAt(2057,iModule); iModule++;
1575 volIDs->AddAt(2058,iModule); iModule++;
1576 volIDs->AddAt(2059,iModule); iModule++;
1577 volIDs->AddAt(2060,iModule); iModule++;
1578 volIDs->AddAt(2061,iModule); iModule++;
1579 volIDs->AddAt(2062,iModule); iModule++;
1580 volIDs->AddAt(2063,iModule); iModule++;
1581 volIDs->AddAt(4112,iModule); iModule++;
1582 volIDs->AddAt(4113,iModule); iModule++;
1583 volIDs->AddAt(4114,iModule); iModule++;
1584 volIDs->AddAt(4115,iModule); iModule++;
1585 volIDs->AddAt(4116,iModule); iModule++;
1586 volIDs->AddAt(4117,iModule); iModule++;
1587 volIDs->AddAt(4118,iModule); iModule++;
1588 volIDs->AddAt(4119,iModule); iModule++;
1589 volIDs->AddAt(4120,iModule); iModule++;
1590 volIDs->AddAt(4121,iModule); iModule++;
1591 volIDs->AddAt(4122,iModule); iModule++;
1592 volIDs->AddAt(4123,iModule); iModule++;
1593 volIDs->AddAt(4124,iModule); iModule++;
1594 volIDs->AddAt(4125,iModule); iModule++;
1595 volIDs->AddAt(4126,iModule); iModule++;
1596 volIDs->AddAt(4127,iModule); iModule++;
1598 if(sectors[2]==1){//--->cSect = 2 <//---
1599 volIDs->AddAt(2064,iModule); iModule++;
1600 volIDs->AddAt(2065,iModule); iModule++;
1601 volIDs->AddAt(2066,iModule); iModule++;
1602 volIDs->AddAt(2067,iModule); iModule++;
1603 volIDs->AddAt(2068,iModule); iModule++;
1604 volIDs->AddAt(2069,iModule); iModule++;
1605 volIDs->AddAt(2070,iModule); iModule++;
1606 volIDs->AddAt(2071,iModule); iModule++;
1607 volIDs->AddAt(4128,iModule); iModule++;
1608 volIDs->AddAt(4129,iModule); iModule++;
1609 volIDs->AddAt(4130,iModule); iModule++;
1610 volIDs->AddAt(4131,iModule); iModule++;
1611 volIDs->AddAt(4132,iModule); iModule++;
1612 volIDs->AddAt(4133,iModule); iModule++;
1613 volIDs->AddAt(4134,iModule); iModule++;
1614 volIDs->AddAt(4135,iModule); iModule++;
1615 volIDs->AddAt(4136,iModule); iModule++;
1616 volIDs->AddAt(4137,iModule); iModule++;
1617 volIDs->AddAt(4138,iModule); iModule++;
1618 volIDs->AddAt(4139,iModule); iModule++;
1619 volIDs->AddAt(4140,iModule); iModule++;
1620 volIDs->AddAt(4141,iModule); iModule++;
1621 volIDs->AddAt(4142,iModule); iModule++;
1622 volIDs->AddAt(4143,iModule); iModule++;
1624 if(sectors[3]==1){//--->cSect = 3 <//---
1625 volIDs->AddAt(2072,iModule); iModule++;
1626 volIDs->AddAt(2073,iModule); iModule++;
1627 volIDs->AddAt(2074,iModule); iModule++;
1628 volIDs->AddAt(2075,iModule); iModule++;
1629 volIDs->AddAt(2076,iModule); iModule++;
1630 volIDs->AddAt(2077,iModule); iModule++;
1631 volIDs->AddAt(2078,iModule); iModule++;
1632 volIDs->AddAt(2079,iModule); iModule++;
1633 volIDs->AddAt(4144,iModule); iModule++;
1634 volIDs->AddAt(4145,iModule); iModule++;
1635 volIDs->AddAt(4146,iModule); iModule++;
1636 volIDs->AddAt(4147,iModule); iModule++;
1637 volIDs->AddAt(4148,iModule); iModule++;
1638 volIDs->AddAt(4149,iModule); iModule++;
1639 volIDs->AddAt(4150,iModule); iModule++;
1640 volIDs->AddAt(4151,iModule); iModule++;
1641 volIDs->AddAt(4152,iModule); iModule++;
1642 volIDs->AddAt(4153,iModule); iModule++;
1643 volIDs->AddAt(4154,iModule); iModule++;
1644 volIDs->AddAt(4155,iModule); iModule++;
1645 volIDs->AddAt(4156,iModule); iModule++;
1646 volIDs->AddAt(4157,iModule); iModule++;
1647 volIDs->AddAt(4158,iModule); iModule++;
1648 volIDs->AddAt(4159,iModule); iModule++;
1650 if(sectors[4]==1){//--->cSect = 4 <//---
1651 volIDs->AddAt(2080,iModule); iModule++;
1652 volIDs->AddAt(2081,iModule); iModule++;
1653 volIDs->AddAt(2082,iModule); iModule++;
1654 volIDs->AddAt(2083,iModule); iModule++;
1655 volIDs->AddAt(2084,iModule); iModule++;
1656 volIDs->AddAt(2085,iModule); iModule++;
1657 volIDs->AddAt(2086,iModule); iModule++;
1658 volIDs->AddAt(2087,iModule); iModule++;
1659 volIDs->AddAt(4160,iModule); iModule++;
1660 volIDs->AddAt(4161,iModule); iModule++;
1661 volIDs->AddAt(4162,iModule); iModule++;
1662 volIDs->AddAt(4163,iModule); iModule++;
1663 volIDs->AddAt(4164,iModule); iModule++;
1664 volIDs->AddAt(4165,iModule); iModule++;
1665 volIDs->AddAt(4166,iModule); iModule++;
1666 volIDs->AddAt(4167,iModule); iModule++;
1667 volIDs->AddAt(4168,iModule); iModule++;
1668 volIDs->AddAt(4169,iModule); iModule++;
1669 volIDs->AddAt(4170,iModule); iModule++;
1670 volIDs->AddAt(4171,iModule); iModule++;
1671 volIDs->AddAt(4172,iModule); iModule++;
1672 volIDs->AddAt(4173,iModule); iModule++;
1673 volIDs->AddAt(4174,iModule); iModule++;
1674 volIDs->AddAt(4175,iModule); iModule++;
1676 if(sectors[5]==1){//--->cSect = 5 <//---
1677 volIDs->AddAt(2088,iModule); iModule++;
1678 volIDs->AddAt(2089,iModule); iModule++;
1679 volIDs->AddAt(2090,iModule); iModule++;
1680 volIDs->AddAt(2091,iModule); iModule++;
1681 volIDs->AddAt(2092,iModule); iModule++;
1682 volIDs->AddAt(2093,iModule); iModule++;
1683 volIDs->AddAt(2094,iModule); iModule++;
1684 volIDs->AddAt(2095,iModule); iModule++;
1685 volIDs->AddAt(4176,iModule); iModule++;
1686 volIDs->AddAt(4177,iModule); iModule++;
1687 volIDs->AddAt(4178,iModule); iModule++;
1688 volIDs->AddAt(4179,iModule); iModule++;
1689 volIDs->AddAt(4180,iModule); iModule++;
1690 volIDs->AddAt(4181,iModule); iModule++;
1691 volIDs->AddAt(4182,iModule); iModule++;
1692 volIDs->AddAt(4183,iModule); iModule++;
1693 volIDs->AddAt(4184,iModule); iModule++;
1694 volIDs->AddAt(4185,iModule); iModule++;
1695 volIDs->AddAt(4186,iModule); iModule++;
1696 volIDs->AddAt(4187,iModule); iModule++;
1697 volIDs->AddAt(4188,iModule); iModule++;
1698 volIDs->AddAt(4189,iModule); iModule++;
1699 volIDs->AddAt(4190,iModule); iModule++;
1700 volIDs->AddAt(4191,iModule); iModule++;
1702 if(sectors[6]==1){//--->cSect = 6 <//---
1703 volIDs->AddAt(2096,iModule); iModule++;
1704 volIDs->AddAt(2097,iModule); iModule++;
1705 volIDs->AddAt(2098,iModule); iModule++;
1706 volIDs->AddAt(2099,iModule); iModule++;
1707 volIDs->AddAt(2100,iModule); iModule++;
1708 volIDs->AddAt(2101,iModule); iModule++;
1709 volIDs->AddAt(2102,iModule); iModule++;
1710 volIDs->AddAt(2103,iModule); iModule++;
1711 volIDs->AddAt(4192,iModule); iModule++;
1712 volIDs->AddAt(4193,iModule); iModule++;
1713 volIDs->AddAt(4194,iModule); iModule++;
1714 volIDs->AddAt(4195,iModule); iModule++;
1715 volIDs->AddAt(4196,iModule); iModule++;
1716 volIDs->AddAt(4197,iModule); iModule++;
1717 volIDs->AddAt(4198,iModule); iModule++;
1718 volIDs->AddAt(4199,iModule); iModule++;
1719 volIDs->AddAt(4200,iModule); iModule++;
1720 volIDs->AddAt(4201,iModule); iModule++;
1721 volIDs->AddAt(4202,iModule); iModule++;
1722 volIDs->AddAt(4203,iModule); iModule++;
1723 volIDs->AddAt(4204,iModule); iModule++;
1724 volIDs->AddAt(4205,iModule); iModule++;
1725 volIDs->AddAt(4206,iModule); iModule++;
1726 volIDs->AddAt(4207,iModule); iModule++;
1728 if(sectors[7]==1){ //--->cSect = 7 <//---
1729 volIDs->AddAt(2104,iModule); iModule++;
1730 volIDs->AddAt(2105,iModule); iModule++;
1731 volIDs->AddAt(2106,iModule); iModule++;
1732 volIDs->AddAt(2107,iModule); iModule++;
1733 volIDs->AddAt(2108,iModule); iModule++;
1734 volIDs->AddAt(2109,iModule); iModule++;
1735 volIDs->AddAt(2110,iModule); iModule++;
1736 volIDs->AddAt(2111,iModule); iModule++;
1737 volIDs->AddAt(4208,iModule); iModule++;
1738 volIDs->AddAt(4209,iModule); iModule++;
1739 volIDs->AddAt(4210,iModule); iModule++;
1740 volIDs->AddAt(4211,iModule); iModule++;
1741 volIDs->AddAt(4212,iModule); iModule++;
1742 volIDs->AddAt(4213,iModule); iModule++;
1743 volIDs->AddAt(4214,iModule); iModule++;
1744 volIDs->AddAt(4215,iModule); iModule++;
1745 volIDs->AddAt(4216,iModule); iModule++;
1746 volIDs->AddAt(4217,iModule); iModule++;
1747 volIDs->AddAt(4218,iModule); iModule++;
1748 volIDs->AddAt(4219,iModule); iModule++;
1749 volIDs->AddAt(4220,iModule); iModule++;
1750 volIDs->AddAt(4221,iModule); iModule++;
1751 volIDs->AddAt(4222,iModule); iModule++;
1752 volIDs->AddAt(4223,iModule); iModule++;
1754 if(sectors[8]==1){//--->cSect = 8 <//---
1755 volIDs->AddAt(2112,iModule); iModule++;
1756 volIDs->AddAt(2113,iModule); iModule++;
1757 volIDs->AddAt(2114,iModule); iModule++;
1758 volIDs->AddAt(2115,iModule); iModule++;
1759 volIDs->AddAt(2116,iModule); iModule++;
1760 volIDs->AddAt(2117,iModule); iModule++;
1761 volIDs->AddAt(2118,iModule); iModule++;
1762 volIDs->AddAt(2119,iModule); iModule++;
1763 volIDs->AddAt(4224,iModule); iModule++;
1764 volIDs->AddAt(4225,iModule); iModule++;
1765 volIDs->AddAt(4226,iModule); iModule++;
1766 volIDs->AddAt(4227,iModule); iModule++;
1767 volIDs->AddAt(4228,iModule); iModule++;
1768 volIDs->AddAt(4229,iModule); iModule++;
1769 volIDs->AddAt(4230,iModule); iModule++;
1770 volIDs->AddAt(4231,iModule); iModule++;
1771 volIDs->AddAt(4232,iModule); iModule++;
1772 volIDs->AddAt(4233,iModule); iModule++;
1773 volIDs->AddAt(4234,iModule); iModule++;
1774 volIDs->AddAt(4235,iModule); iModule++;
1775 volIDs->AddAt(4236,iModule); iModule++;
1776 volIDs->AddAt(4237,iModule); iModule++;
1777 volIDs->AddAt(4238,iModule); iModule++;
1778 volIDs->AddAt(4239,iModule); iModule++;
1780 if(sectors[9]==1){//--->cSect = 9 <//---
1781 volIDs->AddAt(2120,iModule); iModule++;
1782 volIDs->AddAt(2121,iModule); iModule++;
1783 volIDs->AddAt(2122,iModule); iModule++;
1784 volIDs->AddAt(2123,iModule); iModule++;
1785 volIDs->AddAt(2124,iModule); iModule++;
1786 volIDs->AddAt(2125,iModule); iModule++;
1787 volIDs->AddAt(2126,iModule); iModule++;
1788 volIDs->AddAt(2127,iModule); iModule++;
1789 volIDs->AddAt(4240,iModule); iModule++;
1790 volIDs->AddAt(4241,iModule); iModule++;
1791 volIDs->AddAt(4242,iModule); iModule++;
1792 volIDs->AddAt(4243,iModule); iModule++;
1793 volIDs->AddAt(4244,iModule); iModule++;
1794 volIDs->AddAt(4245,iModule); iModule++;
1795 volIDs->AddAt(4246,iModule); iModule++;
1796 volIDs->AddAt(4247,iModule); iModule++;
1797 volIDs->AddAt(4248,iModule); iModule++;
1798 volIDs->AddAt(4249,iModule); iModule++;
1799 volIDs->AddAt(4250,iModule); iModule++;
1800 volIDs->AddAt(4251,iModule); iModule++;
1801 volIDs->AddAt(4252,iModule); iModule++;
1802 volIDs->AddAt(4253,iModule); iModule++;
1803 volIDs->AddAt(4254,iModule); iModule++;
1804 volIDs->AddAt(4255,iModule); iModule++;
1811 //____________________________________________________________________________
1812 TArrayI* AliITSResidualsAnalysis::GetITSLayersVolids(Int_t layers[6]) const
1815 // This method gets the volID Array for the chosen layers.
1816 // You have to pass an array with a 1 for each selected layer.
1817 // i.e. layers[6] = {1,1,0,0,1,1} -> SPD + SSD
1820 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1822 Int_t size=0,last=0;
1824 // evaluates the size of the array
1825 for(Int_t i=0;i<6;i++) if(layers[i]==1) size+=AliGeomManager::LayerSize(i+1);
1828 printf("Error: no layer selected");
1832 TArrayI *volids = new TArrayI(size);
1834 // fills the volId array only for the chosen layers
1835 for(Int_t ilayer=1;ilayer<7;ilayer++){
1837 if(layers[ilayer-1]!=1) continue;
1839 for(Int_t imod=0;imod<AliGeomManager::LayerSize(ilayer);imod++){
1840 volids->AddAt(AliGeomManager::LayerToVolUID(ilayer,imod),last);
1849 //____________________________________________________________________________
1850 void AliITSResidualsAnalysis::GetTrackDirClusterCov(AliTrackPoint *point,Double_t &phi,Double_t &lambda,Double_t &lambda2,Double_t &alpha,Double_t &xovery,Double_t &zovery) const
1857 const Float_t *covvector=point->GetCov();
1858 cov(0,0)=covvector[0];
1859 cov(1,0)=cov(0,1)=covvector[1];
1860 cov(2,0)=cov(0,2)=covvector[2];
1861 cov(1,1)=covvector[3];
1862 cov(1,2)=cov(2,1)=covvector[4];
1863 cov(2,2)=covvector[5];
1865 Double_t determinant=cov.Determinant();
1866 if(determinant!=0.){
1868 TVectorD eigenvalues(3);
1869 const TMatrixDSymEigen keigen(cov);
1870 eigenvalues=keigen.GetEigenValues();
1871 vect=keigen.GetEigenVectors();
1872 Double_t mainvect[3];
1873 mainvect[0]=vect(0,0);
1874 mainvect[1]=vect(1,0);
1875 mainvect[2]=vect(2,0);
1876 if(mainvect[1]!=0.){
1877 xovery=mainvect[0]/mainvect[1];
1878 zovery=mainvect[2]/mainvect[1];
1885 mainvect[0]=-1.*mainvect[0];
1886 mainvect[1]=-1.*mainvect[1];
1887 mainvect[2]=-1.*mainvect[2];
1889 lambda2=TMath::ATan2(TMath::Sqrt(mainvect[0]*mainvect[0]+mainvect[2]*mainvect[2]),mainvect[1])*TMath::RadToDeg();
1890 lambda=TMath::ATan2(mainvect[2],TMath::Sqrt(mainvect[0]*mainvect[0]+mainvect[1]*mainvect[1]))*TMath::RadToDeg();
1891 phi=TMath::ATan2(mainvect[0],mainvect[2])*TMath::RadToDeg();
1892 alpha=TMath::ATan2(mainvect[1],mainvect[0])*TMath::RadToDeg();
1894 else printf("determinant =0!, skip this point \n");
1899 //____________________________________________________________________________
1900 void AliITSResidualsAnalysis::CalculateResiduals(const TArrayI *volids,
1901 const TArrayI *volidsfit,
1902 AliGeomManager::ELayerID layerRangeMin,
1903 AliGeomManager::ELayerID layerRangeMax,
1906 // CalculateResiduals for a set of detector volumes.
1907 // Tracks are fitted only within
1908 // the range defined by the user
1909 // (by layerRangeMin and layerRangeMax)
1910 // or within the set of volidsfit
1911 // Repeat the procedure 'iterations' times
1913 Int_t nVolIds = volids->GetSize();
1914 if (nVolIds == 0) { AliError("Volume IDs array is empty!"); return; }
1916 // Load only the tracks with at least one
1917 // space point in the set of volume (volids)
1920 //AliAlignmentTracks::SetPointsFilename(GetFileNameTrackPoints());
1921 AliAlignmentTracks::BuildIndex();
1923 ListVolUsed(fPointsTree,fArrayIndex,fLastIndex);
1924 AliTrackPointArray **points;
1927 LoadPoints(volids, points,pointsDim);
1929 Int_t nArrays = fPointsTree->GetEntries();
1931 if (nArrays == 0){ AliError("Points array is empty!"); return; }
1932 AliTrackFitter *fitter = CreateFitter();
1937 for (Int_t iArray = 0; iArray < nArrays; iArray++){
1939 cout<<"Investigating "<<iArray<<"/"<<nArrays<<endl;
1941 if (!points[iArray]){
1942 cout<<" Skipping: "<<iArray<<endl;
1946 fitter->SetTrackPointArray(points[iArray],kTRUE); // Watch out, problems
1952 if(fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE){
1954 cout<<"->BAD: "<<iArray<<endl;
1956 } //else cout<<"->GOOD: "<<iArray<<endl;
1958 AliTrackPointArray *pVolId,*pTrack;
1960 fitter->GetTrackResiduals(pVolId,pTrack);
1961 FillResidualsH(pVolId,pTrack);
1965 cout<<" -> nVolIds: "<<nVolIds<<endl;
1966 cout<<" -> Non-Fitted tracks: "<<ecount<<"/"<<totcount<<endl;
1968 UnloadPoints(pointsDim,points);
1969 SaveHists(3,outname);
1977 //______________________________________________________________________________
1978 void AliITSResidualsAnalysis::ProcessVolumes(Int_t fit,
1981 TString misalignmentFile,
1986 // This function process the AliTrackPoints and volID (into residuals)
1989 // setting up geometry and the trackpoints file
1990 if(!gGeoManager) AliGeomManager::LoadGeometry(GetFileNameGeometry());
1992 SetPointsFilename(GetFileNameTrackPoints());
1994 // creating some tools
1995 AliTrackFitter *fitter;
1997 fitter = new AliTrackFitterKalman();
1998 }else fitter = new AliTrackFitterRieman();
2000 fitter->SetMinNPoints(minPoints);
2002 SetTrackFitter(fitter);
2004 if(misalignmentFile=="")printf("NO FAKE MISALIGNMENT\n");
2006 Bool_t misal=Misalign(misalignmentFile,"ITSAlignObjs");
2008 printf("PROBLEM WITH FAKE MISALIGNMENT!");
2013 CalculateResiduals(volIDs,volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,outname);