[u/mrichter/AliRoot.git] / STEER / AliESDV0MI.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $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
//-------------------------------------------------------------------------

#include <Riostream.h>
#include <TMath.h>

#include "AliESDV0MI.h"

ClassImp(AliESDV0MI)

AliESDV0MIParams  AliESDV0MI::fgkParams;


AliESDV0MI::AliESDV0MI() :
  AliESDv0(),
  fParamP(),
  fParamM(),
  fID(0),
  fDist1(-1),
  fDist2(-1),
  fRr(-1),
  fStatus(0),
  fRow0(-1),
  fDistNorm(0),
  fDistSigma(0),
  fChi2Before(0),
  fNBefore(0),
  fChi2After(0),
  fNAfter(0),
  fPointAngleFi(0),
  fPointAngleTh(0),
  fPointAngle(0)
{
  //
  //Dafault constructor
  //
  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 (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
}


Double_t AliESDV0MI::GetEffectiveSigmaAP0(){
  //
  // effective Sigma parameterization of point angle resolution 
  //
  Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
			     fParamM.GetParameter()[4]*fParamM.GetParameter()[4]);
  Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
  sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
  sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
  sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
  return sigmaAPE;
}


Double_t  AliESDV0MI::GetMinimaxSigmaAP0(){
  //
  // calculate mini-max effective sigma of point angle resolution
  //
  //compv0->fTree->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)
{
  //
  // set probabilities
  //
  fCausality[0] = pb0;     // probability - track 0 exist before vertex
  fCausality[1] = pb1;     // probability - track 1 exist before vertex
  fCausality[2] = pa0;     // probability - track 0 exist close after vertex
  fCausality[3] = pa1;     // probability - track 1 exist close after vertex
}
void  AliESDV0MI::SetClusters(Int_t *clp, Int_t *clm)
{
  //
  // Set its clusters indexes
  //
  for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i]; 
  for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i]; 
}


void AliESDV0MI::SetP(const AliExternalTrackParam & paramp)  {
  //
  // set track +
  //
  fParamP   = paramp;
}

void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
  //
  //set track -
  //
  fParamM = paramm;
}
  
void AliESDV0MI::SetRp(const Double_t *rp){
  //
  // set pid +
  //
  for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
}

void AliESDV0MI::SetRm(const Double_t *rm){
  //
  // set pid -
  //
  for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
}


void  AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
{
  //
  // set PID hypothesy
  //
  // norm PID to 1
  Float_t sump =0;
  Float_t summ =0;
  for (Int_t i=0;i<5;i++){
    fRP[i]=pidp[i];
    sump+=fRP[i];
    fRM[i]=pidm[i];
    summ+=fRM[i];
  }
  for (Int_t i=0;i<5;i++){
    fRP[i]/=sump;
    fRM[i]/=summ;
  }
}

Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
  //
  //
  //
  //
  return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
}

Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
  //
  // calculate effective mass
  //
  const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
		      4.93599999999999983e-01, 9.38270000000000048e-01};
  if (p1>4) return -1;
  if (p2>4) return -1;
  Float_t mass1 = kpmass[p1]; 
  Float_t mass2 = kpmass[p2];   
  Double_t *m1 = fPP;
  Double_t *m2 = fPM;
  //
  //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
  //  m1 = fPM;
  //  m2 = fPP;
  //}
  //
  Float_t e1    = TMath::Sqrt(mass1*mass1+
                              m1[0]*m1[0]+
                              m1[1]*m1[1]+
                              m1[2]*m1[2]);
  Float_t e2    = TMath::Sqrt(mass2*mass2+
                              m2[0]*m2[0]+
                              m2[1]*m2[1]+
                              m2[2]*m2[2]);  
  Float_t mass =  
    (m2[0]+m1[0])*(m2[0]+m1[0])+
    (m2[1]+m1[1])*(m2[1]+m1[1])+
    (m2[2]+m1[2])*(m2[2]+m1[2]);
  
  mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
  return mass;
}
Commit	Line	Data
51ad6848	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	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	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	//-------------------------------------------------------------------------
c1e38247	19	//
	20	// Implementation of the ESD V0MI vertex class
	21	// This class is part of the Event Data Summary
	22	// set of classes and contains information about
	23	// V0 kind vertexes generated by a neutral particle
	24	// Numerical part - AliHelix functionality used
	25	//
	26	// Likelihoods for Point angle, DCA and Causality defined => can be used as cut parameters
	27	// HIGHLY recomended
	28	//
	29	// Quality information can be used as additional cut variables
	30	//
51ad6848	31	// Origin: Marian Ivanov marian.ivanov@cern.ch
	32	//-------------------------------------------------------------------------
	33
	34	#include <Riostream.h>
	35	#include <TMath.h>
0703142d	36
51ad6848	37	#include "AliESDV0MI.h"
51ad6848	38
	39	ClassImp(AliESDV0MI)
	40
c1e38247	41	AliESDV0MIParams AliESDV0MI::fgkParams;
	42
	43
90e48c0c	44	AliESDV0MI::AliESDV0MI() :
	45	AliESDv0(),
	46	fParamP(),
	47	fParamM(),
	48	fID(0),
	49	fDist1(-1),
	50	fDist2(-1),
	51	fRr(-1),
	52	fStatus(0),
	53	fRow0(-1),
	54	fDistNorm(0),
	55	fDistSigma(0),
	56	fChi2Before(0),
	57	fNBefore(0),
	58	fChi2After(0),
	59	fNAfter(0),
	60	fPointAngleFi(0),
	61	fPointAngleTh(0),
	62	fPointAngle(0)
	63	{
51ad6848	64	//
	65	//Dafault constructor
	66	//
eaacfdf5	67	for (Int_t i=0;i<5;i++){
	68	fRP[i]=fRM[i]=0;
	69	}
	70	fLab[0]=fLab[1]=-1;
	71	fIndex[0]=fIndex[1]=-1;
6605de26	72	for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
eaacfdf5	73	fNormDCAPrim[0]=fNormDCAPrim[1]=0;
	74	for (Int_t i=0;i<3;i++){fPP[i]=fPM[i]=fXr[i]=fAngle[i]=0;}
	75	for (Int_t i=0;i<3;i++){fOrder[i]=0;}
	76	for (Int_t i=0;i<4;i++){fCausality[i]=0;}
81e97e0d	77	}
81e97e0d	78
c1e38247	79	Double_t AliESDV0MI::GetSigmaY(){
	80	//
	81	// return sigmay in y at vertex position using covariance matrix
	82	//
	83	const Double_t * cp = fParamP.GetCovariance();
	84	const Double_t * cm = fParamM.GetCovariance();
c9ec41e8	85	Double_t sigmay = cp[0]+cm[0]+ cp[5](fParamP.GetX()-fRr)(fParamP.GetX()-fRr)+ cm[5](fParamM.GetX()-fRr)(fParamM.GetX()-fRr);
c1e38247	86	return (sigmay>0) ? TMath::Sqrt(sigmay):100;
	87	}
	88
	89	Double_t AliESDV0MI::GetSigmaZ(){
	90	//
	91	// return sigmay in y at vertex position using covariance matrix
	92	//
	93	const Double_t * cp = fParamP.GetCovariance();
	94	const Double_t * cm = fParamM.GetCovariance();
c9ec41e8	95	Double_t sigmaz = cp[2]+cm[2]+ cp[9](fParamP.GetX()-fRr)(fParamP.GetX()-fRr)+ cm[9](fParamM.GetX()-fRr)(fParamM.GetX()-fRr);
c1e38247	96	return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
	97	}
	98
	99	Double_t AliESDV0MI::GetSigmaD0(){
	100	//
	101	// Sigma parameterization using covariance matrix
	102	//
	103	// sigma of distance between two tracks in vertex position
	104	// sigma of DCA is proportianal to sigmaD0
	105	// factor 2 difference is explained by the fact that the DCA is calculated at the position
	106	// where the tracks as closest together ( not exact position of the vertex)
	107	//
	108	const Double_t * cp = fParamP.GetCovariance();
	109	const Double_t * cm = fParamM.GetCovariance();
	110	Double_t sigmaD0 = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
c9ec41e8	111	sigmaD0 += ((fParamP.GetX()-fRr)(fParamP.GetX()-fRr))(cp[5]+cp[9]);
c9ec41e8	112	sigmaD0 += ((fParamM.GetX()-fRr)(fParamM.GetX()-fRr))(cm[5]+cm[9]);
c1e38247	113	return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
	114	}
	115
	116
	117	Double_t AliESDV0MI::GetSigmaAP0(){
	118	//
	119	//Sigma parameterization using covariance matrices
	120	//
	121	Double_t prec = TMath::Sqrt((fPM[0]+fPP[0])*(fPM[0]+fPP[0])
	122	+(fPM[1]+fPP[1])*(fPM[1]+fPP[1])
	123	+(fPM[2]+fPP[2])*(fPM[2]+fPP[2]));
	124	Double_t normp = TMath::Sqrt(fPP[0]fPP[0]+fPP[1]fPP[1]+fPP[2]*fPP[2])/prec; // fraction of the momenta
	125	Double_t normm = TMath::Sqrt(fPM[0]fPM[0]+fPM[1]fPM[1]+fPM[2]*fPM[2])/prec;
	126	const Double_t * cp = fParamP.GetCovariance();
	127	const Double_t * cm = fParamM.GetCovariance();
	128	Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0; // minimal part
	129	sigmaAP0 += (cp[5]+cp[9])(normpnormp)+(cm[5]+cm[9])(normmnormm); // angular resolution part
	130	Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01); // vertex position part
	131	sigmaAP0 += 0.5sigmaAP1sigmaAP1;
	132	return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
	133	}
	134
	135	Double_t AliESDV0MI::GetEffectiveSigmaD0(){
	136	//
	137	// minimax - effective Sigma parameterization
	138	// p12 effective curvature and v0 radius postion used as parameters
	139	//
	140	Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
	141	fParamM.GetParameter()[4]*fParamM.GetParameter()[4]);
	142	Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)fgkParams.fPSigmaCoefDEp12*p12;
	143	sigmaED0*= sigmaED0;
	144	sigmaED0*= sigmaED0;
	145	sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
	146	return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
	147	}
	148
	149
	150	Double_t AliESDV0MI::GetEffectiveSigmaAP0(){
	151	//
	152	// effective Sigma parameterization of point angle resolution
	153	//
	154	Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
	155	fParamM.GetParameter()[4]*fParamM.GetParameter()[4]);
	156	Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
	157	sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
	158	sigmaAPE= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APEp12);
	159	sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
	160	return sigmaAPE;
	161	}
	162
	163
	164	Double_t AliESDV0MI::GetMinimaxSigmaAP0(){
	165	//
	166	// calculate mini-max effective sigma of point angle resolution
	167	//
	168	//compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)0.5,1.5SigmaAPE0),0.5*SigmaAPE0+0.003)");
	169	Double_t effectiveSigma = GetEffectiveSigmaAP0();
	170	Double_t sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
	171	sigmaMMAP = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
	172	sigmaMMAP = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
	173	return sigmaMMAP;
	174	}
	175	Double_t AliESDV0MI::GetMinimaxSigmaD0(){
	176	//
177	// calculate mini-max sigma of dca resolution
178	//
179	//compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)0.5,1.5SigmaDE0),0.5*SigmaDE0)");
180	Double_t effectiveSigma = GetEffectiveSigmaD0();
181	Double_t sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
182	sigmaMMD0 = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
183	sigmaMMD0 = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
184	return sigmaMMD0;
185	}
186
187
188	Double_t AliESDV0MI::GetLikelihoodAP(Int_t mode0, Int_t mode1){
189	//
190	// get likelihood for point angle
191	//
192	Double_t sigmaAP = 0.007; //default sigma
193	switch (mode0){
194	case 0:
195	sigmaAP = GetSigmaAP0(); // mode 0 - covariance matrix estimates used
196	break;
197	case 1:
198	sigmaAP = GetEffectiveSigmaAP0(); // mode 1 - effective sigma used
199	break;
200	case 2:
201	sigmaAP = GetMinimaxSigmaAP0(); // mode 2 - minimax sigma
202	break;
203	}
204	Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);
205	//normalized point angle, restricted - because of overflow problems in Exp
206	Double_t likelihood = 0;
207	switch(mode1){
208	case 0:
209	likelihood = TMath::Exp(-0.5apNormapNorm);
210	// one component
211	break;
212	case 1:
213	likelihood = (TMath::Exp(-0.5apNormapNorm)+0.5* TMath::Exp(-0.25apNormapNorm))/1.5;
214	// two components
215	break;
216	case 2:
217	likelihood = (TMath::Exp(-0.5apNormapNorm)+0.5* TMath::Exp(-0.25apNormapNorm)+0.25TMath::Exp(-0.125apNorm*apNorm))/1.75;
218	// three components
219	break;
220	}
221	return likelihood;
222	}
223
224	Double_t AliESDV0MI::GetLikelihoodD(Int_t mode0, Int_t mode1){
225	//
226	// get likelihood for DCA
227	//
228	Double_t sigmaD = 0.03; //default sigma
229	switch (mode0){
230	case 0:
231	sigmaD = GetSigmaD0(); // mode 0 - covariance matrix estimates used
232	break;
233	case 1:
234	sigmaD = GetEffectiveSigmaD0(); // mode 1 - effective sigma used
235	break;
236	case 2:
237	sigmaD = GetMinimaxSigmaD0(); // mode 2 - minimax sigma
238	break;
239	}
240	Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
241	//normalized point angle, restricted - because of overflow problems in Exp
242	Double_t likelihood = 0;
243	switch(mode1){
244	case 0:
245	likelihood = TMath::Exp(-2.*dNorm);
246	// one component
247	break;
248	case 1:
249	likelihood = (TMath::Exp(-2.dNorm)+0.5 TMath::Exp(-dNorm))/1.5;
250	// two components
251	break;
252	case 2:
253	likelihood = (TMath::Exp(-2.dNorm)+0.5 TMath::Exp(-dNorm)+0.25TMath::Exp(-0.5dNorm))/1.75;
254	// three components
255	break;
256	}
257	return likelihood;
258
259	}
260
261	Double_t AliESDV0MI::GetLikelihoodC(Int_t mode0, Int_t /mode1/){
262	//
263	// get likelihood for Causality
264	// !!! Causality variables defined in AliITStrackerMI !!!
265	// when more information was available
266	//
267	Double_t likelihood = 0.5;
268	Double_t minCausal = TMath::Min(fCausality[0],fCausality[1]);
269	Double_t maxCausal = TMath::Max(fCausality[0],fCausality[1]);
270	// minCausal = TMath::Max(minCausal,0.5*maxCausal);
271	//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");
272
273	switch(mode0){
274	case 0:
275	//normalization
276	likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
277	break;
278	case 1:
279	likelihood = TMath::Power(1.05-(2(0.8-TMath::Exp(-maxCausal))+(2(0.8-TMath::Exp(-minCausal))))*0.5,4.);
280	break;
281	}
282	return likelihood;
283
284	}
81e97e0d	285
	286	void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
	287	{
	288	//
	289	// set probabilities
	290	//
	291	fCausality[0] = pb0; // probability - track 0 exist before vertex
	292	fCausality[1] = pb1; // probability - track 1 exist before vertex
	293	fCausality[2] = pa0; // probability - track 0 exist close after vertex
	294	fCausality[3] = pa1; // probability - track 1 exist close after vertex
51ad6848	295	}
6605de26	296	void AliESDV0MI::SetClusters(Int_t clp, Int_t clm)
	297	{
	298	//
	299	// Set its clusters indexes
	300	//
	301	for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i];
	302	for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i];
	303	}
	304
51ad6848	305
	306	void AliESDV0MI::SetP(const AliExternalTrackParam & paramp) {
	307	//
81e97e0d	308	// set track +
51ad6848	309	//
	310	fParamP = paramp;
	311	}
	312
	313	void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
	314	//
81e97e0d	315	//set track -
51ad6848	316	//
51ad6848	317	fParamM = paramm;
51ad6848	318	}
51ad6848	319
81e97e0d	320	void AliESDV0MI::SetRp(const Double_t *rp){
	321	//
	322	// set pid +
	323	//
	324	for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
	325	}
	326
	327	void AliESDV0MI::SetRm(const Double_t *rm){
	328	//
	329	// set pid -
	330	//
	331	for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
	332	}
	333
	334
51ad6848	335	void AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
	336	{
	337	//
	338	// set PID hypothesy
	339	//
	340	// norm PID to 1
	341	Float_t sump =0;
	342	Float_t summ =0;
	343	for (Int_t i=0;i<5;i++){
	344	fRP[i]=pidp[i];
	345	sump+=fRP[i];
	346	fRM[i]=pidm[i];
	347	summ+=fRM[i];
	348	}
	349	for (Int_t i=0;i<5;i++){
	350	fRP[i]/=sump;
	351	fRM[i]/=summ;
	352	}
	353	}
	354
	355	Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
	356	//
	357	//
	358	//
	359	//
	360	return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
	361	}
	362
	363	Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
	364	//
	365	// calculate effective mass
	366	//
0703142d	367	const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
51ad6848	368	4.93599999999999983e-01, 9.38270000000000048e-01};
	369	if (p1>4) return -1;
	370	if (p2>4) return -1;
0703142d	371	Float_t mass1 = kpmass[p1];
0703142d	372	Float_t mass2 = kpmass[p2];
51ad6848	373	Double_t *m1 = fPP;
	374	Double_t *m2 = fPM;
	375	//
6605de26	376	//if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
	377	// m1 = fPM;
	378	// m2 = fPP;
	379	//}
51ad6848	380	//
	381	Float_t e1 = TMath::Sqrt(mass1*mass1+
	382	m1[0]*m1[0]+
	383	m1[1]*m1[1]+
	384	m1[2]*m1[2]);
	385	Float_t e2 = TMath::Sqrt(mass2*mass2+
	386	m2[0]*m2[0]+
	387	m2[1]*m2[1]+
	388	m2[2]*m2[2]);
	389	Float_t mass =
	390	(m2[0]+m1[0])*(m2[0]+m1[0])+
	391	(m2[1]+m1[1])*(m2[1]+m1[1])+
	392	(m2[2]+m1[2])*(m2[2]+m1[2]);
	393
	394	mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
	395	return mass;
	396	}
	397