#include <TChain.h>
#include <TFile.h>
+#include <TVector3.h>
+#include <TSystem.h>
#include "AliAlignmentTracks.h"
#include "AliTrackPointArray.h"
-#include "AliAlignObjAngles.h"
+#include "AliAlignObjParams.h"
#include "AliTrackFitterRieman.h"
#include "AliTrackResidualsChi2.h"
-#include "AliESD.h"
+#include "AliESDEvent.h"
#include "AliLog.h"
ClassImp(AliAlignmentTracks)
fMisalignObjs(0),
fTrackFitter(0),
fMinimizer(0),
- fDoUpdate(kTRUE)
+ fDoUpdate(kTRUE),
+ fCovIsUsed(kFALSE)
{
// Default constructor
InitIndex();
fMisalignObjs(0),
fTrackFitter(0),
fMinimizer(0),
- fDoUpdate(kTRUE)
+ fDoUpdate(kTRUE),
+ fCovIsUsed(kFALSE)
{
// Constructor in the case
// the user provides an already
fMisalignObjs(0),
fTrackFitter(0),
fMinimizer(0),
- fDoUpdate(kTRUE)
+ fDoUpdate(kTRUE),
+ fCovIsUsed(kFALSE)
{
// Constructor in the case
// the user provides a single ESD file
}
}
-//________________________________________________________________________
+//________________________________________________________________________
void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
Int_t minITSpts,
Bool_t cuts,
+ Float_t minAngleWrtITSModulePlanes,
Float_t minMom,Float_t maxMom,
Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
Float_t minSinTheta,Float_t maxSinTheta)
if (!fESDChain) return;
- AliESD *esd = 0;
- fESDChain->SetBranchAddress("ESD",&esd);
+ AliESDEvent *esd = new AliESDEvent();
+ esd->ReadFromTree(fESDChain);
AliESDfriend *esdf = 0;
fESDChain->SetBranchStatus("ESDfriend*",1);
fESDChain->SetBranchAddress("ESDfriend.",&esdf);
// Open the output file
- if (fPointsFilename.Data() == "") {
+ if (fPointsFilename.IsNull()) {
AliWarning("Incorrect output filename!");
return;
}
Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
}
- // UInt_t status = AliESDtrack::kITSpid;
- // status|=AliESDtrack::kTPCpid;
- // status|=AliESDtrack::kTRDpid;
- // if ((track->GetStatus() & status) != status) continue;
AliTrackPoint point;
array = track->GetTrackPointArray();
if(onlyITS) {
- array2 = new AliTrackPointArray(track->GetNcls(0));
+ Bool_t layerOK[6]={kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
Int_t ipt,volId,modId,layerId;
Int_t jpt=0;
for(ipt=0; ipt<array->GetNPoints(); ipt++) {
array->GetPoint(point,ipt);
volId = point.GetVolumeID();
layerId = AliGeomManager::VolUIDToLayer(volId,modId);
- if(layerId<7) {
- array2->AddPoint(jpt,&point);
- jpt++;
+ if(layerId>6) continue;
+ // check minAngleWrtITSModulePlanes
+ if(cuts) {
+ Double_t p[3]; track->GetDirection(p);
+ TVector3 pvec(p[0],p[1],p[2]);
+ Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
+ TVector3 normvec(rot[1],rot[4],rot[7]);
+ Double_t angle = pvec.Angle(normvec);
+ if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
+ angle = 0.5*TMath::Pi()-angle;
+ if(angle<minAngleWrtITSModulePlanes) {
+ layerOK[layerId-1]=kFALSE;
+ continue;
+ }
}
+ jpt++;
+ }
+ if(jpt < minITSpts) continue;
+ array2 = new AliTrackPointArray(jpt);
+ jpt=0;
+ for(ipt=0; ipt<array->GetNPoints(); ipt++) {
+ array->GetPoint(point,ipt);
+ volId = point.GetVolumeID();
+ layerId = AliGeomManager::VolUIDToLayer(volId,modId);
+ if(layerId>6 || !layerOK[layerId-1]) continue;
+ array2->AddPoint(jpt,&point);
+ jpt++;
}
} // end if(onlyITS)
}
pointsFile->Close();
+
+ return;
}
-//______________________________________________________________________________
+//_____________________________________________________________________________
void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
- Int_t minITSpts,
- Bool_t cuts,
+ Int_t minITSpts,Float_t maxMatchingAngle,
+ Bool_t cuts,
+ Float_t minAngleWrtITSModulePlanes,
Float_t minMom,Float_t maxMom,
Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
Float_t minSinTheta,Float_t maxSinTheta)
if (!fESDChain) return;
- AliESD *esd = 0;
- fESDChain->SetBranchAddress("ESD",&esd);
+ AliESDEvent *esd = new AliESDEvent();
+ esd->ReadFromTree(fESDChain);
AliESDfriend *esdf = 0;
fESDChain->SetBranchStatus("ESDfriend*",1);
fESDChain->SetBranchAddress("ESDfriend.",&esdf);
// Open the output file
- if (fPointsFilename.Data() == "") {
+ if (fPointsFilename.IsNull()) {
AliWarning("Incorrect output filename!");
return;
}
esd->SetESDfriend(esdf); //Attach the friend to the ESD
Int_t ntracks = esd->GetNumberOfTracks();
+ if(ntracks<2) continue;
+ Int_t *goodtracksArray = new Int_t[ntracks];
+ Float_t *phiArray = new Float_t[ntracks];
+ Float_t *thetaArray = new Float_t[ntracks];
Int_t ngt=0;
- Int_t goodtracks[10];
for (Int_t itrack=0; itrack < ntracks; itrack++) {
AliESDtrack * track = esd->GetTrack(itrack);
if (!track) continue;
if(track->GetNcls(0) < minITSpts) continue;
+ Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
+ Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
if(cuts) {
if(track->GetP()<minMom || track->GetP()>maxMom) continue;
- Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
+ Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
- Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
+ Float_t sintheta = TMath::Sin(theta);
if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
}
- if(ngt<10) goodtracks[ngt]=itrack;
+ goodtracksArray[ngt]=itrack;
+ phiArray[ngt]=phi;
+ thetaArray[ngt]=theta;
ngt++;
}
- if(ngt!=2) continue; // this can be changed to check that the two tracks match
-
- AliESDtrack * track1 = esd->GetTrack(goodtracks[0]);
- AliESDtrack * track2 = esd->GetTrack(goodtracks[1]);
+ if(ngt<2) {
+ delete [] goodtracksArray; goodtracksArray=0;
+ delete [] phiArray; phiArray=0;
+ delete [] thetaArray; thetaArray=0;
+ continue;
+ }
- Int_t ntotpts;
- if(onlyITS) {
- ntotpts = track1->GetNcls(0)+track2->GetNcls(0);
- } else {
- ntotpts = track1->GetTrackPointArray()->GetNPoints()+
- track2->GetTrackPointArray()->GetNPoints();
+ // check matching of the two tracks from the muon
+ Float_t min = 10000000.;
+ Int_t good1 = -1, good2 = -1;
+ for(Int_t itr1=0; itr1<ngt-1; itr1++) {
+ for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
+ Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
+ if(deltatheta>maxMatchingAngle) continue;
+ Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
+ if(deltaphi>maxMatchingAngle) continue;
+ //printf("%f %f %f %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
+ if(deltatheta+deltaphi<min) {
+ min=deltatheta+deltaphi;
+ good1 = goodtracksArray[itr1];
+ good2 = goodtracksArray[itr2];
+ }
+ }
}
- array2 = new AliTrackPointArray(ntotpts);
+
+ delete [] goodtracksArray; goodtracksArray=0;
+ delete [] phiArray; phiArray=0;
+ delete [] thetaArray; thetaArray=0;
+
+ if(good1<0) continue;
+
+ AliESDtrack * track1 = esd->GetTrack(good1);
+ AliESDtrack * track2 = esd->GetTrack(good2);
+
AliTrackPoint point;
Int_t ipt,volId,modId,layerId;
Int_t jpt=0;
+ Bool_t layerOK[6][2];
+ for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kTRUE;
array = track1->GetTrackPointArray();
for(ipt=0; ipt<array->GetNPoints(); ipt++) {
array->GetPoint(point,ipt);
volId = point.GetVolumeID();
layerId = AliGeomManager::VolUIDToLayer(volId,modId);
if(layerId>6) continue;
+ // check minAngleWrtITSModulePlanes
+ if(cuts) {
+ Double_t p[3]; track1->GetDirection(p);
+ TVector3 pvec(p[0],p[1],p[2]);
+ Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
+ TVector3 normvec(rot[1],rot[4],rot[7]);
+ Double_t angle = pvec.Angle(normvec);
+ if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
+ angle = 0.5*TMath::Pi()-angle;
+ if(angle<minAngleWrtITSModulePlanes) {
+ layerOK[layerId-1][0]=kFALSE;
+ continue;
+ }
+ }
}
- array2->AddPoint(jpt,&point);
jpt++;
}
array = track2->GetTrackPointArray();
volId = point.GetVolumeID();
layerId = AliGeomManager::VolUIDToLayer(volId,modId);
if(layerId>6) continue;
+ // check minAngleWrtITSModulePlanes
+ if(cuts) {
+ Double_t p[3]; track2->GetDirection(p);
+ TVector3 pvec(p[0],p[1],p[2]);
+ Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
+ TVector3 normvec(rot[1],rot[4],rot[7]);
+ Double_t angle = pvec.Angle(normvec);
+ if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
+ angle = 0.5*TMath::Pi()-angle;
+ if(angle<minAngleWrtITSModulePlanes) {
+ layerOK[layerId-1][0]=kFALSE;
+ continue;
+ }
+ }
+ }
+ jpt++;
+ }
+
+ if(jpt < 2*minITSpts) continue;
+ array2 = new AliTrackPointArray(jpt);
+ jpt=0;
+ array = track1->GetTrackPointArray();
+ for(ipt=0; ipt<array->GetNPoints(); ipt++) {
+ array->GetPoint(point,ipt);
+ if(onlyITS) {
+ volId = point.GetVolumeID();
+ layerId = AliGeomManager::VolUIDToLayer(volId,modId);
+ if(layerId>6 || !layerOK[layerId-1][0]) continue;
+ }
+ array2->AddPoint(jpt,&point);
+ jpt++;
+ }
+ array = track2->GetTrackPointArray();
+ for(ipt=0; ipt<array->GetNPoints(); ipt++) {
+ array->GetPoint(point,ipt);
+ if(onlyITS) {
+ volId = point.GetVolumeID();
+ layerId = AliGeomManager::VolUIDToLayer(volId,modId);
+ if(layerId>6 || !layerOK[layerId-1][1]) continue;
}
array2->AddPoint(jpt,&point);
jpt++;
}
pointsFile->Close();
+ return;
}
//______________________________________________________________________________
// Dummy object is created in order
// to initialize the volume paths
- AliAlignObjAngles alobj;
+ AliAlignObjParams alobj;
fPointsFile = TFile::Open(fPointsFilename);
if (!fPointsFile || !fPointsFile->IsOpen()) {
//______________________________________________________________________________
Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
{
- // Read alignment object from a file
+ // Read alignment object from a file: update the alignobj already present with the one in the file
// To be replaced by a call to CDB
- AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
+
+ if(gSystem->AccessPathName(alignObjFileName,kFileExists)){
+ printf("Wrong AlignObjs File Name \n");
+ return kFALSE;
+ }
- return kFALSE;
+ TFile *fRealign=TFile::Open(alignObjFileName);
+ if (!fRealign || !fRealign->IsOpen()) {
+ AliError(Form("Could not open Align Obj File file %s !",alignObjFileName));
+ return kFALSE;
+ }
+ printf("Getting TClonesArray \n");
+ TClonesArray *clnarray=(TClonesArray*)fRealign->Get(arrayName);
+ Int_t size=clnarray->GetSize();
+ UShort_t volid;
+
+ for(Int_t ivol=0;ivol<size;ivol++){
+ AliAlignObjParams *a=(AliAlignObjParams*)clnarray->At(ivol);
+ volid=a->GetVolUID();
+ Int_t iModule;
+ AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
+ if(iLayer<AliGeomManager::kFirstLayer||iLayer>AliGeomManager::kSSD2)continue;
+ printf("Updating volume: %d ,layer: %d module: %d \n",volid,iLayer,iModule);
+ *fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule] *= *a;
+ }
+
+ delete clnarray;
+ fRealign->Close();
+ return kTRUE;
}
//______________________________________________________________________________
fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
- fAlignObjs[iLayer][iModule] = new AliAlignObjAngles(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
+ fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
}
}
}
fAlignObjs = 0;
}
-void AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
- AliGeomManager::ELayerID lastLayer,
- AliGeomManager::ELayerID layerRangeMin,
- AliGeomManager::ELayerID layerRangeMax,
- Int_t iterations)
+Bool_t AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
+ AliGeomManager::ELayerID lastLayer,
+ AliGeomManager::ELayerID layerRangeMin,
+ AliGeomManager::ELayerID layerRangeMax,
+ Int_t iterations)
{
// Align detector volumes within
// a given layer range
}
}
+ Bool_t result = kFALSE;
while (iterations > 0) {
- AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
+ if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
iterations--;
}
+ return result;
}
//______________________________________________________________________________
-void AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
- AliGeomManager::ELayerID layerRangeMin,
- AliGeomManager::ELayerID layerRangeMax,
- Int_t iterations)
+Bool_t AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
+ AliGeomManager::ELayerID layerRangeMin,
+ AliGeomManager::ELayerID layerRangeMax,
+ Int_t iterations)
{
// Align detector volumes within
// a given layer.
volIds.AddAt(volId,iModule);
}
+ Bool_t result = kFALSE;
while (iterations > 0) {
- AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
+ if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
iterations--;
}
+ return result;
}
//______________________________________________________________________________
-void AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
+Bool_t AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
Int_t iterations)
{
// Align single detector volume to
TArrayI volIdsFit(1);
volIdsFit.AddAt(volIdFit,0);
+ Bool_t result = kFALSE;
while (iterations > 0) {
- AlignVolumes(&volIds,&volIdsFit);
+ if (!(result = AlignVolumes(&volIds,&volIdsFit))) break;
iterations--;
}
+ return result;
}
//______________________________________________________________________________
-void AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
+Bool_t AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
AliGeomManager::ELayerID layerRangeMin,
AliGeomManager::ELayerID layerRangeMax,
Int_t iterations)
Int_t nVolIds = volids->GetSize();
if (nVolIds == 0) {
AliError("Volume IDs array is empty!");
- return;
+ return kFALSE;
}
-
+
// Load only the tracks with at least one
// space point in the set of volume (volids)
BuildIndex();
AliTrackPointArray **points;
+ Int_t pointsdim;
// Start the iterations
+ Bool_t result = kFALSE;
while (iterations > 0) {
- Int_t nArrays = LoadPoints(volids, points);
- if (nArrays == 0) return;
+ Int_t nArrays = LoadPoints(volids, points,pointsdim);
+ if (nArrays == 0) return kFALSE;
AliTrackResiduals *minimizer = CreateMinimizer();
minimizer->SetNTracks(nArrays);
fitter->GetTrackResiduals(pVolId,pTrack);
minimizer->AddTrackPointArrays(pVolId,pTrack);
}
- minimizer->Minimize();
+ if (!(result = minimizer->Minimize())) break;
// Update the alignment object(s)
if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
*alignObj *= *minimizer->GetAlignObj();
- alignObj->Print("");
+ if(iterations==1)alignObj->Print("");
}
- UnloadPoints(nArrays, points);
+ UnloadPoints(pointsdim, points);
iterations--;
}
+ return result;
}
//______________________________________________________________________________
-Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points)
+Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points,Int_t &pointsdim)
{
// Load track point arrays with at least
// one space point in a given set of detector
if (nArrays == 0) {
AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
- points = 0;
+ points = 0x0;
return 0;
}
fPointsTree->SetBranchAddress("SP", &array);
// Allocate the pointer to the space-point arrays
+ pointsdim=nArrays;
points = new AliTrackPointArray*[nArrays];
for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
// Get tree entry
Int_t entry = (*index)[iArrayId];
- if (indexUsed[entry] == kTRUE) continue;
+ if (indexUsed[entry] == kTRUE) {
+ nArrays--;
+ continue;
+ }
fPointsTree->GetEvent(entry);
if (!array) {
AliWarning("Wrong space point array index!");
}
// End of misalignment
+
AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
- alignObj->Transform(p);
+ UShort_t volp=p.GetVolumeID();
+ Bool_t found=kFALSE;
+ if(fCovIsUsed){
+ for (Int_t iVol = 0; iVol < nVolIds; iVol++) {
+ UShort_t vol = (*volids)[iVol];
+ if(volp==vol){
+ alignObj->Transform(p,kFALSE);
+ found=kTRUE;
+ break;
+ }
+ }
+ }
+ if(!found)alignObj->Transform(p,fCovIsUsed);
points[iArray]->AddPoint(iPoint,&p);
}
iArray++;
{
AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
alObj->GetVolUID(layerId,modId);
+ if(layerId<AliGeomManager::kFirstLayer) {
+ AliWarning(Form("Alignment object is ignored: %s",alObj->GetSymName()));
+ continue;
+ }
fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
}
return kTRUE;
}
+
+
+//________________________________________________
+void AliAlignmentTracks::WriteRealignObjArray(TString outfilename,AliGeomManager::ELayerID layerRangeMin,AliGeomManager::ELayerID layerRangeMax){
+
+ Int_t last=0;
+ TClonesArray *clonesArray=new TClonesArray("AliAlignObjParams",2200);
+ TClonesArray &alo=*clonesArray;
+ for (Int_t iLayer = layerRangeMin-AliGeomManager::kFirstLayer; iLayer <= (layerRangeMax - AliGeomManager::kFirstLayer);iLayer++) {
+
+ for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
+
+ AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
+ new(alo[last])AliAlignObjParams(*alignObj);
+ last++;
+ }
+ }
+ TFile *file=new TFile(outfilename.Data(),"RECREATE");
+ file->cd();
+
+ alo.Write("ITSAlignObjs",TObject::kSingleKey);
+ file->Close();
+
+
+ delete clonesArray;
+ return;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff