#include "TDatabasePDG.h"
#include "TMath.h"
#include "AliLog.h"
-#include "TClonesArray.h"
+#include "TObjArray.h"
+#include "TDecompLU.h"
#include "AliMFTAnalysisTools.h"
Bool_t AliMFTAnalysisTools::ExtrapAODMuonToZ(AliAODTrack *muon, Double_t z, Double_t xy[2], TLorentzVector &kinem, TMatrixD &cov) {
+ // Extrapolate muon to a given z providing the corresponding (x,y) position and updating kinematics and covariance matrix
+
if (!(muon->Pz()!=0)) return kFALSE;
AliMUONTrackParam *param = new AliMUONTrackParam();
//====================================================================================================================================================
-Double_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
+Bool_t AliMFTAnalysisTools::ExtrapAODMuonToXY(AliAODTrack *muon, Double_t xy[2], Double_t &zFinal, TLorentzVector &kinem, TMatrixD &cov) {
+
+ // Find the point of closest approach between the muon and the direction parallel to the z-axis defined by the given (x,y)
+ // Provide the z of the above point as weel as the updated kinematics and covariance matrix
+
+ // We look for the above-defined PCA
+
+ AliMUONTrackParam *param = new AliMUONTrackParam();
+ param -> SetNonBendingCoor(muon->XAtDCA());
+ param -> SetBendingCoor(muon->YAtDCA());
+ param -> SetZ(AliMFTConstants::fZEvalKinem);
+ param -> SetNonBendingSlope(muon->Px()/muon->Pz());
+ param -> SetBendingSlope(muon->Py()/muon->Pz());
+ param -> SetInverseBendingMomentum( muon->Charge() * (1./muon->Pz()) / (TMath::Sqrt(1+TMath::Power(muon->Py()/muon->Pz(),2))) );
+
+ // here we want to understand in which direction we have to search the minimum...
+
+ Double_t step = 1.; // initial step, in cm
+ Double_t startPoint = 0.;
+
+ Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
+
+ TVector3 **points = new TVector3*[2]; // points[0] for the muon, points[1] for the direction parallel to the z-axis defined by the given (x,y)
+
+ for (Int_t i=0; i<3; i++) {
+ AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
+ points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
+ points[1] = new TVector3(xy[0],xy[1],z[i]);
+ r[i] = GetDistanceBetweenPoints(points,2);
+ for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
+ }
+
+ Int_t researchDirection = 0;
+
+ if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
+ else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
+ else if (r[0]<r[1] && r[1]>r[2]) {
+ printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima)\n");
+ delete param;
+ delete points;
+ return kFALSE;
+ }
+
+ while (TMath::Abs(researchDirection)>0.5) {
+
+ if (researchDirection>0) {
+ z[0] = z[1];
+ z[1] = z[2];
+ z[2] = z[1]+researchDirection*step;
+ }
+ else {
+ z[2] = z[1];
+ z[1] = z[0];
+ z[0] = z[1]+researchDirection*step;
+ }
+ if (TMath::Abs(z[0])>900.) {
+ printf("E-AliMFTAnalysisTools::ExtrapAODMuonToXY: Point of closest approach cannot be found (no minima in the fiducial region)\n");
+ delete param;
+ delete points;
+ return kFALSE;
+ }
+
+ for (Int_t i=0; i<3; i++) {
+ AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
+ points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
+ points[1] = new TVector3(xy[0],xy[1],z[i]);
+ r[i] = GetDistanceBetweenPoints(points,2);
+ for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
+ }
+
+ researchDirection=0;
+ if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
+ else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
+
+ }
+
+ // now we now that the minimum is between z[0] and z[2] and we search for it
+
+ step *= 0.5;
+ while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
+ z[0] = z[1]-step;
+ z[2] = z[1]+step;
+ for (Int_t i=0; i<3; i++) {
+ AliMUONTrackExtrap::ExtrapToZ(param, z[i]);
+ points[0] = new TVector3(param->GetNonBendingCoor(),param->GetBendingCoor(),z[i]);
+ points[1] = new TVector3(xy[0],xy[1],z[i]);
+ r[i] = GetDistanceBetweenPoints(points,2);
+ for (Int_t iMu=0; iMu<2; iMu++) delete points[iMu];
+ }
+ if (r[0]<r[1]) z[1] = z[0];
+ else if (r[2]<r[1]) z[1] = z[2];
+ else step *= 0.5;
+ }
+
+ zFinal = z[1];
+
+ Double_t xyMuon[2] = {0};
+ ExtrapAODMuonToZ(muon, zFinal, xyMuon, kinem, cov);
+
+ return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::GetAODMuonOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
Double_t xy[2] = {0};
ExtrapAODMuonToZ(muon, zv, xy);
- return TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
+ offset = TMath::Sqrt((xv-xy[0])*(xv-xy[0]) + (yv-xy[1])*(yv-xy[1]));
+
+ return kTRUE;
}
//====================================================================================================================================================
-Double_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv) {
+Bool_t AliMFTAnalysisTools::GetAODMuonOffsetZ(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
+
+ Double_t xy[2] = {xv, yv};
+ Double_t zFinal = 0;
+ TLorentzVector kinem(0,0,0,0);
+ TMatrixD cov(5,5);
+
+ ExtrapAODMuonToXY(muon, xy, zFinal, kinem, cov);
+
+ offset = TMath::Abs(zFinal - zv);
+
+ return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double_t xv, Double_t yv, Double_t zv, Double_t &offset) {
Double_t xy[2] = {0};
TLorentzVector kinem(0,0,0,0);
covCoordinates(1,0) = cov(2,0);
covCoordinates(1,1) = cov(2,2);
- TMatrixD covCoordinatesInverse = covCoordinates.Invert();
+ if (TDecompLU::InvertLU(covCoordinates,covCoordinates.GetTol(),0)) {
- Double_t dX = xy[0] - xv;
- Double_t dY = xy[1] - yv;
-
- Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
- dY*dY*covCoordinatesInverse(1,1) +
- 2.*dX*dY*covCoordinatesInverse(0,1)));
+ TMatrixD covCoordinatesInverse = covCoordinates;
+ Double_t dX = xy[0] - xv;
+ Double_t dY = xy[1] - yv;
+
+ offset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
+ dY*dY*covCoordinatesInverse(1,1) +
+ 2.*dX*dY*covCoordinatesInverse(0,1)));
+
+ return kTRUE;
- return weightedOffset;
+ }
+
+ return kFALSE;
}
//====================================================================================================================================================
-Bool_t AliMFTAnalysisTools::CalculatePCA(TClonesArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
+Bool_t AliMFTAnalysisTools::CalculatePCA(AliAODDimuon *dimuon, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
+
+ TObjArray *muons = new TObjArray();
+ muons -> Add(dimuon->GetMu(0));
+ muons -> Add(dimuon->GetMu(1));
+
+ Bool_t result = CalculatePCA(muons, pca, pcaQuality, kinem);
+ delete muons;
+ return result;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::CalculatePCA(TObjArray *muons, Double_t *pca, Double_t &pcaQuality, TLorentzVector &kinem) {
const Int_t nMuons = muons->GetEntriesFast();
if (nMuons<2 || nMuons>AliMFTConstants::fNMaxMuonsForPCA) {
}
Double_t fXPointOfClosestApproach=0, fYPointOfClosestApproach=0, fZPointOfClosestApproach=0;
-
- AliAODTrack *muon[nMuons];
- AliMUONTrackParam *param[nMuons];
+
+ AliAODTrack *muon[AliMFTConstants::fNMaxMuonsForPCA] = {0};
+ AliMUONTrackParam *param[AliMFTConstants::fNMaxMuonsForPCA] = {0};
// Finding AliMUONTrackParam objects for each muon
for (Int_t iMu=0; iMu<nMuons; iMu++) {
muon[iMu] = (AliAODTrack*) muons->At(iMu);
- if (TMath::Abs(muon[iMu]->Pz())<1.e-6) return kFALSE;
+ if (TMath::Abs(muon[iMu]->Pz())<1.e-6) {
+ for(Int_t i=0;i<iMu;i++) delete param[i];
+ return kFALSE;
+ }
param[iMu] = new AliMUONTrackParam();
param[iMu] -> SetNonBendingCoor(muon[iMu]->XAtDCA());
param[iMu] -> SetBendingCoor(muon[iMu]->YAtDCA());
- param[iMu] -> SetZ(0.);
+ param[iMu] -> SetZ(AliMFTConstants::fZEvalKinem);
param[iMu] -> SetNonBendingSlope(muon[iMu]->Px()/muon[iMu]->Pz());
param[iMu] -> SetBendingSlope(muon[iMu]->Py()/muon[iMu]->Pz());
param[iMu] -> SetInverseBendingMomentum( muon[iMu]->Charge() * (1./muon[iMu]->Pz()) / (TMath::Sqrt(1+TMath::Power(muon[iMu]->Py()/muon[iMu]->Pz(),2))) );
}
-
+
// here we want to understand in which direction we have to search the minimum...
Double_t step = 1.; // initial step, in cm
Double_t startPoint = 0.;
-
+
Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
- TVector3 **points = new TVector3*[nMuons];
+ TVector3 **points = new TVector3*[AliMFTConstants::fNMaxMuonsForPCA];
for (Int_t i=0; i<3; i++) {
for (Int_t iMu=0; iMu<nMuons; iMu++) {
AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
- }
+ }
r[i] = GetDistanceBetweenPoints(points,nMuons);
+ for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu];
}
-
+
Int_t researchDirection = 0;
if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
else if (r[0]<r[1] && r[1]>r[2]) {
printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima)\n");
+ for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
+ delete points;
return kFALSE;
- }
-
+ }
+
while (TMath::Abs(researchDirection)>0.5) {
-
+
if (researchDirection>0) {
z[0] = z[1];
z[1] = z[2];
}
if (TMath::Abs(z[0])>900.) {
printf("E-AliMFTAnalysisTools::CalculatePCA: Point of closest approach cannot be found for dimuon (no minima in the fiducial region)\n");
+ for (Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
+ delete points;
return kFALSE;
}
- for (Int_t iMu=0; iMu<nMuons; iMu++) delete points[iMu];
-
for (Int_t i=0; i<3; i++) {
for (Int_t iMu=0; iMu<nMuons; iMu++) {
- AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
- points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
+ AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
+ points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
}
r[i] = GetDistanceBetweenPoints(points,nMuons);
+ for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
}
researchDirection=0;
if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1; // towards z positive
else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1; // towards z negative
-
+
}
-
+
// now we now that the minimum is between z[0] and z[2] and we search for it
step *= 0.5;
z[2] = z[1]+step;
for (Int_t i=0; i<3; i++) {
for (Int_t iMu=0; iMu<nMuons; iMu++) {
- AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
- points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
+ AliMUONTrackExtrap::ExtrapToZ(param[iMu], z[i]);
+ points[iMu] = new TVector3(param[iMu]->GetNonBendingCoor(),param[iMu]->GetBendingCoor(),z[i]);
}
r[i] = GetDistanceBetweenPoints(points,nMuons);
+ for (Int_t iMu=0;iMu<nMuons;iMu++) delete points[iMu];
}
if (r[0]<r[1]) z[1] = z[0];
else if (r[2]<r[1]) z[1] = z[2];
else step *= 0.5;
}
- delete [] points;
-
// Once z of minimum is found, we evaluate the x and y coordinates by averaging over the contributing tracks
fZPointOfClosestApproach = z[1];
pca[2] = fZPointOfClosestApproach;
// Evaluating the kinematics of the N-muon
-
+
Double_t pTot[3] = {0};
Double_t ene = 0.;
Double_t massMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
}
kinem.SetPxPyPzE(pTot[0], pTot[1], pTot[2], ene);
-
+
// Evaluating the PCA quality of the N-muon
-
+
Double_t sum=0.,squareSum=0.;
for (Int_t iMu=0; iMu<nMuons; iMu++) {
- Double_t wOffset = AliMFTAnalysisTools::GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
+ Double_t wOffset = 0;
+ if (!GetAODMuonWeightedOffset(muon[iMu],fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach, wOffset)) {
+ for(Int_t jMu=0;jMu<nMuons;jMu++) delete param[jMu];
+ delete points;
+ return kFALSE;
+ }
Double_t f = TMath::Exp(-0.5 * wOffset);
sum += f;
squareSum += f*f;
}
if (sum > 0.) pcaQuality = (sum-squareSum/sum) / (nMuons-1);
else pcaQuality = 0.;
-
+
+ for(Int_t iMu=0;iMu<nMuons;iMu++) delete param[iMu];
+ delete points;
return kTRUE;
}
//====================================================================================================================================================
+Bool_t AliMFTAnalysisTools::IsCorrectMatch(AliAODTrack *muon) {
+
+ for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) if (IsWrongCluster(muon, iPlane)) return kFALSE;
+ return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+TString AliMFTAnalysisTools::GetGenerator(Int_t label, AliAODMCHeader* header) {
+
+ // get the name of the generator that produced a given particle
+
+ Int_t partCounter = 0;
+ TList *genHeaders = header->GetCocktailHeaders();
+ Int_t nGenHeaders = genHeaders->GetEntries();
+
+ for (Int_t i=0; i<nGenHeaders; i++){
+ AliGenEventHeader *gh = (AliGenEventHeader*) genHeaders->At(i);
+ TString genName = gh->GetName();
+ Int_t nPart = gh->NProduced();
+ if (label>=partCounter && label<(partCounter+nPart)) return genName;
+ partCounter += nPart;
+ }
+
+ TString empty="";
+ return empty;
+
+}
+
+//====================================================================================================================================================
+
+void AliMFTAnalysisTools::GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen) {
+
+ // method to check if a track comes from a given generator
+
+ Int_t label = TMath::Abs(track->GetLabel());
+ nameGen = GetGenerator(label,header);
+
+ // In case the particle is not primary nameGen will contain blank spaces. In this case, we search backward for the primary which originated the chain
+
+ while (nameGen.IsWhitespace()) {
+ AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(label);
+ if (!mcPart) {
+ printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: No valid AliAODMCParticle at label %i\n",label);
+ break;
+ }
+ Int_t motherLabel = mcPart->GetMother();
+ if (motherLabel < 0) {
+ printf("AliMFTAnalysisTools::GetTrackPrimaryGenerator - BREAK: Reached primary particle without valid mother\n");
+ break;
+ }
+ label = motherLabel;
+ nameGen = GetGenerator(label,header);
+ }
+
+ return;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTAnalysisTools::IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC){
+
+ // method to check if a track comes from the signal event or from the underlying Hijing event
+
+ TString nameGen;
+
+ GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
+
+ if (nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
+
+ return kTRUE;
+}
+
+//====================================================================================================================================================
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff