#include <TChain.h>
#include <TFile.h>
+#include <TVector3.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)
}
}
-//______________________________________________________________________________
-void AliAlignmentTracks::ProcessESD()
+
+//________________________________________________________________________
+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)
{
// Analyzes and filters ESD tracks
// Stores the selected track space points
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);
TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
const AliTrackPointArray *array = 0;
- pointsTree->Branch("SP","AliTrackPointArray", &array);
+ AliTrackPointArray *array2 = 0;
+ if(onlyITS) { // only ITS AliTrackPoints
+ pointsTree->Branch("SP","AliTrackPointArray", &array2);
+ } else {
+ pointsTree->Branch("SP","AliTrackPointArray", &array);
+ }
+
Int_t ievent = 0;
while (fESDChain->GetEntry(ievent++)) {
if (!esd) break;
for (Int_t itrack=0; itrack < ntracks; itrack++) {
AliESDtrack * track = esd->GetTrack(itrack);
if (!track) continue;
-
- // UInt_t status = AliESDtrack::kITSpid;
-// status|=AliESDtrack::kTPCpid;
-// status|=AliESDtrack::kTRDpid;
-// if ((track->GetStatus() & status) != status) continue;
- if (track->GetP() < 0.5) continue;
+ if(track->GetNcls(0) < minITSpts) continue;
+ if(cuts) {
+ if(track->GetP()<minMom || track->GetP()>maxMom) continue;
+ Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
+ if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
+ Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
+ if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
+ }
+ AliTrackPoint point;
array = track->GetTrackPointArray();
+
+ if(onlyITS) {
+ 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>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[0],0.);
+ 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)
+
pointsTree->Fill();
}
}
}
pointsFile->Close();
+
+ return;
+}
+
+//_____________________________________________________________________________
+void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
+ 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)
+{
+ // Analyzes and filters ESD tracks
+ // Merges inward and outward tracks in one single track
+ // Stores the selected track space points
+ // into the output file
+
+ if (!fESDChain) return;
+
+ 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() == "") {
+ AliWarning("Incorrect output filename!");
+ return;
+ }
+
+ TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
+ if (!pointsFile || !pointsFile->IsOpen()) {
+ AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
+ return;
+ }
+
+ TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
+ const AliTrackPointArray *array = 0;
+ AliTrackPointArray *array2 = 0;
+ pointsTree->Branch("SP","AliTrackPointArray", &array2);
+
+ Int_t ievent = 0;
+ while (fESDChain->GetEntry(ievent++)) {
+ if (!esd) break;
+
+ 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;
+ 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(phi));
+ if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
+ Float_t sintheta = TMath::Sin(theta);
+ if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
+ }
+ goodtracksArray[ngt]=itrack;
+ phiArray[ngt]=phi;
+ thetaArray[ngt]=theta;
+ ngt++;
+ }
+
+ if(ngt<2) {
+ delete [] goodtracksArray; goodtracksArray=0;
+ delete [] phiArray; phiArray=0;
+ delete [] thetaArray; thetaArray=0;
+ continue;
+ }
+
+ // 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];
+ }
+ }
+ }
+
+ 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);
+ if(onlyITS) {
+ 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[0],0.);
+ 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++;
+ }
+ 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) 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[0],0.);
+ 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++;
+ }
+
+ pointsTree->Fill();
+ }
+
+ if (!pointsTree->Write()) {
+ AliWarning("Can't write the tree with track point arrays!");
+ return;
+ }
+
+ 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()) {
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);
}
}
}
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);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff