X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliESDV0MI.cxx;h=4d9ae028ec860b67fb9a86cb02ae5f9393be820e;hb=f753671bb4e4862a3cbdd49fada4e6def4e90de9;hp=4d6d9c5e0679ffac7b1ae5e01d3c875b6e009fb8;hpb=6605de26edef346181dbf28e61bad6d5d151c0ee;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliESDV0MI.cxx b/STEER/AliESDV0MI.cxx index 4d6d9c5e067..4d9ae028ec8 100644 --- a/STEER/AliESDV0MI.cxx +++ b/STEER/AliESDV0MI.cxx @@ -16,6 +16,18 @@ /* $Id$ */ //------------------------------------------------------------------------- +// +// Implementation of the ESD V0MI vertex class +// This class is part of the Event Data Summary +// set of classes and contains information about +// V0 kind vertexes generated by a neutral particle +// Numerical part - AliHelix functionality used +// +// Likelihoods for Point angle, DCA and Causality defined => can be used as cut parameters +// HIGHLY recomended +// +// Quality information can be used as additional cut variables +// // Origin: Marian Ivanov marian.ivanov@cern.ch //------------------------------------------------------------------------- @@ -28,6 +40,9 @@ ClassImp(AliESDV0MI) +AliESDV0MIParams AliESDV0MI::fgkParams; + + AliESDV0MI::AliESDV0MI() : AliESDv0(), fParamP(), @@ -51,11 +66,225 @@ AliESDV0MI::AliESDV0MI() : // //Dafault constructor // - for (Int_t i=0;i<4;i++){fCausality[i]=0;} + for (Int_t i=0;i<5;i++){ + fRP[i]=fRM[i]=0; + } + fLab[0]=fLab[1]=-1; + fIndex[0]=fIndex[1]=-1; for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;} + fNormDCAPrim[0]=fNormDCAPrim[1]=0; + for (Int_t i=0;i<3;i++){fPP[i]=fPM[i]=fXr[i]=fAngle[i]=0;} + for (Int_t i=0;i<3;i++){fOrder[i]=0;} + for (Int_t i=0;i<4;i++){fCausality[i]=0;} +} + +Double_t AliESDV0MI::GetSigmaY(){ + // + // return sigmay in y at vertex position using covariance matrix + // + const Double_t * cp = fParamP.GetCovariance(); + const Double_t * cm = fParamM.GetCovariance(); + Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamM.GetX()-fRr)*(fParamM.GetX()-fRr); + return (sigmay>0) ? TMath::Sqrt(sigmay):100; +} + +Double_t AliESDV0MI::GetSigmaZ(){ + // + // return sigmay in y at vertex position using covariance matrix + // + const Double_t * cp = fParamP.GetCovariance(); + const Double_t * cm = fParamM.GetCovariance(); + Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamM.GetX()-fRr)*(fParamM.GetX()-fRr); + return (sigmaz>0) ? TMath::Sqrt(sigmaz):100; +} + +Double_t AliESDV0MI::GetSigmaD0(){ + // + // Sigma parameterization using covariance matrix + // + // sigma of distance between two tracks in vertex position + // sigma of DCA is proportianal to sigmaD0 + // factor 2 difference is explained by the fact that the DCA is calculated at the position + // where the tracks as closest together ( not exact position of the vertex) + // + const Double_t * cp = fParamP.GetCovariance(); + const Double_t * cm = fParamM.GetCovariance(); + Double_t sigmaD0 = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0; + sigmaD0 += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]); + sigmaD0 += ((fParamM.GetX()-fRr)*(fParamM.GetX()-fRr))*(cm[5]+cm[9]); + return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100; +} + + +Double_t AliESDV0MI::GetSigmaAP0(){ + // + //Sigma parameterization using covariance matrices + // + Double_t prec = TMath::Sqrt((fPM[0]+fPP[0])*(fPM[0]+fPP[0]) + +(fPM[1]+fPP[1])*(fPM[1]+fPP[1]) + +(fPM[2]+fPP[2])*(fPM[2]+fPP[2])); + Double_t normp = TMath::Sqrt(fPP[0]*fPP[0]+fPP[1]*fPP[1]+fPP[2]*fPP[2])/prec; // fraction of the momenta + Double_t normm = TMath::Sqrt(fPM[0]*fPM[0]+fPM[1]*fPM[1]+fPM[2]*fPM[2])/prec; + const Double_t * cp = fParamP.GetCovariance(); + const Double_t * cm = fParamM.GetCovariance(); + Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0; // minimal part + sigmaAP0 += (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm); // angular resolution part + Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01); // vertex position part + sigmaAP0 += 0.5*sigmaAP1*sigmaAP1; + return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100; +} + +Double_t AliESDV0MI::GetEffectiveSigmaD0(){ + // + // minimax - effective Sigma parameterization + // p12 effective curvature and v0 radius postion used as parameters + // + Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+ + fParamM.GetParameter()[4]*fParamM.GetParameter()[4]); + Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12; + sigmaED0*= sigmaED0; + sigmaED0*= sigmaED0; + sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE); + return (sigmaED0fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)"); + Double_t effectiveSigma = GetEffectiveSigmaAP0(); + Double_t sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma); + sigmaMMAP = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma); + sigmaMMAP = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0); + return sigmaMMAP; +} +Double_t AliESDV0MI::GetMinimaxSigmaD0(){ + // + // calculate mini-max sigma of dca resolution + // + //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)"); + Double_t effectiveSigma = GetEffectiveSigmaD0(); + Double_t sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma); + sigmaMMD0 = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma); + sigmaMMD0 = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0); + return sigmaMMD0; +} + + +Double_t AliESDV0MI::GetLikelihoodAP(Int_t mode0, Int_t mode1){ + // + // get likelihood for point angle + // + Double_t sigmaAP = 0.007; //default sigma + switch (mode0){ + case 0: + sigmaAP = GetSigmaAP0(); // mode 0 - covariance matrix estimates used + break; + case 1: + sigmaAP = GetEffectiveSigmaAP0(); // mode 1 - effective sigma used + break; + case 2: + sigmaAP = GetMinimaxSigmaAP0(); // mode 2 - minimax sigma + break; + } + Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.); + //normalized point angle, restricted - because of overflow problems in Exp + Double_t likelihood = 0; + switch(mode1){ + case 0: + likelihood = TMath::Exp(-0.5*apNorm*apNorm); + // one component + break; + case 1: + likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5; + // two components + break; + case 2: + likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm)+0.25*TMath::Exp(-0.125*apNorm*apNorm))/1.75; + // three components + break; + } + return likelihood; +} + +Double_t AliESDV0MI::GetLikelihoodD(Int_t mode0, Int_t mode1){ + // + // get likelihood for DCA + // + Double_t sigmaD = 0.03; //default sigma + switch (mode0){ + case 0: + sigmaD = GetSigmaD0(); // mode 0 - covariance matrix estimates used + break; + case 1: + sigmaD = GetEffectiveSigmaD0(); // mode 1 - effective sigma used + break; + case 2: + sigmaD = GetMinimaxSigmaD0(); // mode 2 - minimax sigma + break; + } + Double_t dNorm = TMath::Min(fDist2/sigmaD,50.); + //normalized point angle, restricted - because of overflow problems in Exp + Double_t likelihood = 0; + switch(mode1){ + case 0: + likelihood = TMath::Exp(-2.*dNorm); + // one component + break; + case 1: + likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5; + // two components + break; + case 2: + likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75; + // three components + break; + } + return likelihood; + +} + +Double_t AliESDV0MI::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/){ + // + // get likelihood for Causality + // !!! Causality variables defined in AliITStrackerMI !!! + // when more information was available + // + Double_t likelihood = 0.5; + Double_t minCausal = TMath::Min(fCausality[0],fCausality[1]); + Double_t maxCausal = TMath::Max(fCausality[0],fCausality[1]); + // minCausal = TMath::Max(minCausal,0.5*maxCausal); + //compv0->fTree->SetAlias("LCausal","(1.05-(2*(0.8-exp(-max(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1])))+2*(0.8-exp(-min(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1]))))/2)**4"); + + switch(mode0){ + case 0: + //normalization + likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.); + break; + case 1: + likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.); + break; + } + return likelihood; + +} + void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1) { // @@ -204,10 +433,26 @@ void AliESDV0MI::Update(Float_t vertex[3]) points = phelix.GetRPHIintersections(mhelix, phase, radius); delta1=10000,delta2=10000; Double_t d1=1000.,d2=10000.; + Double_t err[3],angles[3]; if (points<=0) return; if (points>0){ phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + if (TMath::Abs(fParamP.GetX()-TMath::Sqrt(radius[0])<3) && TMath::Abs(fParamM.GetX()-TMath::Sqrt(radius[0])<3)){ + // if we are close to vertex use error parama + // + err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05 + +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); + err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05 + +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); + + phelix.GetAngle(phase[0][0],mhelix,phase[0][1],angles); + Double_t tfi = TMath::Abs(TMath::Tan(angles[0])); + Double_t tlam = TMath::Abs(TMath::Tan(angles[1])); + err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam)); + err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0]; + phelix.ParabolicDCA2(mhelix,phase[0][0],phase[0][1],radius[0],delta1,err); + } Double_t xd[3],xm[3]; phelix.Evaluate(phase[0][0],xd); mhelix.Evaluate(phase[0][1],xm); @@ -216,6 +461,21 @@ void AliESDV0MI::Update(Float_t vertex[3]) if (points==2){ phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + if (TMath::Abs(fParamP.GetX()-TMath::Sqrt(radius[1])<3) && TMath::Abs(fParamM.GetX()-TMath::Sqrt(radius[1])<3)){ + // if we are close to vertex use error paramatrization + // + err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05 + +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); + err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05 + +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); + + phelix.GetAngle(phase[1][0],mhelix,phase[1][1],angles); + Double_t tfi = TMath::Abs(TMath::Tan(angles[0])); + Double_t tlam = TMath::Abs(TMath::Tan(angles[1])); + err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam)); + err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0]; + phelix.ParabolicDCA2(mhelix,phase[1][0],phase[1][1],radius[1],delta2,err); + } Double_t xd[3],xm[3]; phelix.Evaluate(phase[1][0],xd); mhelix.Evaluate(phase[1][1],xm); @@ -231,6 +491,12 @@ void AliESDV0MI::Update(Float_t vertex[3]) fXr[0] = 0.5*(xd[0]+xm[0]); fXr[1] = 0.5*(xd[1]+xm[1]); fXr[2] = 0.5*(xd[2]+xm[2]); + + Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); + Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); + fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] ); + fXr[1] = (1.-wy)*xd[1]+ wy*xm[1]; + fXr[2] = (1.-wz)*xd[2]+ wz*xm[2]; // phelix.GetMomentum(phase[0][0],fPP); mhelix.GetMomentum(phase[0][1],fPM); @@ -244,6 +510,11 @@ void AliESDV0MI::Update(Float_t vertex[3]) fXr[0] = 0.5*(xd[0]+xm[0]); fXr[1] = 0.5*(xd[1]+xm[1]); fXr[2] = 0.5*(xd[2]+xm[2]); + Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); + Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); + fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] ); + fXr[1] = (1.-wy)*xd[1]+ wy*xm[1]; + fXr[2] = (1.-wz)*xd[2]+ wz*xm[2]; // phelix.GetMomentum(phase[1][0], fPP); mhelix.GetMomentum(phase[1][1], fPM); @@ -257,7 +528,8 @@ void AliESDV0MI::Update(Float_t vertex[3]) Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]}; Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]}; Double_t vnorm2 = v[0]*v[0]+v[1]*v[1]; - Double_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2); + if (TMath::Abs(v[2])>100000) return; + Double_t vnorm3 = TMath::Sqrt(TMath::Abs(v[2]*v[2]+vnorm2)); vnorm2 = TMath::Sqrt(vnorm2); Double_t pnorm2 = p[0]*p[0]+p[1]*p[1]; Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);