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

//-------------------------------------------------------------------------
//    Origin: Marian Ivanov marian.ivanov@cern.ch
//-------------------------------------------------------------------------

#include <Riostream.h>
#include <TMath.h>
#include <TPDGCode.h>
#include "AliESDV0MI.h"
#include "AliHelix.h"


ClassImp(AliESDV0MI)

AliESDV0MI::AliESDV0MI(){
  //
  //Dafault constructor
  //
  fID =0;
  fDist1  =-1;
  fDist2  =-1;
  fRr     =-1;
  fStatus = 0;
  fRow0   =-1;  
}

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

void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
  //
  //set daughter
  //
  fParamM = paramm;

}
  
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 pmass[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 = pmass[p1]; 
  Float_t mass2 = pmass[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;
}

void  AliESDV0MI::Update(Float_t vertex[3])
{
  //
  // updates Kink Info
  //
  Float_t distance1,distance2;
  //
  AliHelix phelix(fParamP);
  AliHelix mhelix(fParamM);    
  //
  //find intersection linear
  //
  Double_t phase[2][2],radius[2];
  Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
  Double_t delta1=10000,delta2=10000;  

  if (points<=0) return;
  if (points>0){
    phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
  }
  if (points==2){    
    phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
  }
  distance1 = TMath::Min(delta1,delta2);
  //
  //find intersection parabolic
  //
  points = phelix.GetRPHIintersections(mhelix, phase, radius);
  delta1=10000,delta2=10000;  
  Double_t d1=1000.,d2=10000.;
  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);
    Double_t xd[3],xm[3];
    phelix.Evaluate(phase[0][0],xd);
    mhelix.Evaluate(phase[0][1],xm);
    d1 = (xd[0]-xm[0])*(xd[0]-xm[0])+(xd[1]-xm[1])*(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]);
  }
  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);
    Double_t xd[3],xm[3];
    phelix.Evaluate(phase[1][0],xd);
    mhelix.Evaluate(phase[1][1],xm);
    d2 = (xd[0]-xm[0])*(xd[0]-xm[0])+(xd[1]-xm[1])*(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]);
  }
  //
  distance2 = TMath::Min(delta1,delta2);
  if (delta1<delta2){
    //get V0 info
    Double_t xd[3],xm[3];
    phelix.Evaluate(phase[0][0],xd);
    mhelix.Evaluate(phase[0][1], xm);
    fXr[0] = 0.5*(xd[0]+xm[0]);
    fXr[1] = 0.5*(xd[1]+xm[1]);
    fXr[2] = 0.5*(xd[2]+xm[2]);
    //
    phelix.GetMomentum(phase[0][0],fPP);
    mhelix.GetMomentum(phase[0][1],fPM);
    phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
    fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
  }
  else{
    Double_t xd[3],xm[3];
    phelix.Evaluate(phase[1][0],xd);
    mhelix.Evaluate(phase[1][1], xm);
    fXr[0] = 0.5*(xd[0]+xm[0]);
    fXr[1] = 0.5*(xd[1]+xm[1]);
    fXr[2] = 0.5*(xd[2]+xm[2]);
    //
    phelix.GetMomentum(phase[1][0], fPP);
    mhelix.GetMomentum(phase[1][1], fPM);
    phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
    fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
  }
  fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
  fDist2 = TMath::Sqrt(distance2);      
  //            
  //
  Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
  Float_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
  vnorm2 = TMath::Sqrt(vnorm2);
  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
  pnorm2 = TMath::Sqrt(pnorm2);  
  fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
  fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
  fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
  //
}
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	//-------------------------------------------------------------------------
	19	// Origin: Marian Ivanov marian.ivanov@cern.ch
	20	//-------------------------------------------------------------------------
	21
	22	#include <Riostream.h>
	23	#include <TMath.h>
	24	#include <TPDGCode.h>
	25	#include "AliESDV0MI.h"
	26	#include "AliHelix.h"
	27
	28
	29	ClassImp(AliESDV0MI)
	30
	31	AliESDV0MI::AliESDV0MI(){
	32	//
	33	//Dafault constructor
	34	//
	35	fID =0;
	36	fDist1 =-1;
	37	fDist2 =-1;
	38	fRr =-1;
	39	fStatus = 0;
	40	fRow0 =-1;
	41	}
	42
	43	void AliESDV0MI::SetP(const AliExternalTrackParam & paramp) {
	44	//
	45	// set mother
	46	//
	47	fParamP = paramp;
	48	}
	49
	50	void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
	51	//
	52	//set daughter
	53	//
	54	fParamM = paramm;
	55
	56	}
	57
	58	void AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
	59	{
	60	//
	61	// set PID hypothesy
	62	//
	63	// norm PID to 1
	64	Float_t sump =0;
65	Float_t summ =0;
66	for (Int_t i=0;i<5;i++){
67	fRP[i]=pidp[i];
68	sump+=fRP[i];
69	fRM[i]=pidm[i];
70	summ+=fRM[i];
71	}
72	for (Int_t i=0;i<5;i++){
73	fRP[i]/=sump;
74	fRM[i]/=summ;
75	}
76	}
77
78	Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
79	//
80	//
81	//
82	//
83	return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
84	}
85
86	Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
87	//
88	// calculate effective mass
89	//
90	const Float_t pmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
91	4.93599999999999983e-01, 9.38270000000000048e-01};
92	if (p1>4) return -1;
93	if (p2>4) return -1;
94	Float_t mass1 = pmass[p1];
95	Float_t mass2 = pmass[p2];
96	Double_t *m1 = fPP;
97	Double_t *m2 = fPM;
98	//
99	if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
100	m1 = fPM;
101	m2 = fPP;
102	}
103	//
104	Float_t e1 = TMath::Sqrt(mass1*mass1+
105	m1[0]*m1[0]+
106	m1[1]*m1[1]+
107	m1[2]*m1[2]);
108	Float_t e2 = TMath::Sqrt(mass2*mass2+
109	m2[0]*m2[0]+
110	m2[1]*m2[1]+
111	m2[2]*m2[2]);
112	Float_t mass =
113	(m2[0]+m1[0])*(m2[0]+m1[0])+
114	(m2[1]+m1[1])*(m2[1]+m1[1])+
115	(m2[2]+m1[2])*(m2[2]+m1[2]);
116
117	mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
118	return mass;
119	}
120
121	void AliESDV0MI::Update(Float_t vertex[3])
122	{
123	//
124	// updates Kink Info
125	//
126	Float_t distance1,distance2;
127	//
128	AliHelix phelix(fParamP);
129	AliHelix mhelix(fParamM);
130	//
131	//find intersection linear
132	//
133	Double_t phase[2][2],radius[2];
134	Int_t points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
135	Double_t delta1=10000,delta2=10000;
b07a036b	136
b07a036b	137	if (points<=0) return;
51ad6848	138	if (points>0){
	139	phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	140	phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	141	phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	142	}
	143	if (points==2){
	144	phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	145	phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	146	phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	147	}
	148	distance1 = TMath::Min(delta1,delta2);
	149	//
	150	//find intersection parabolic
	151	//
	152	points = phelix.GetRPHIintersections(mhelix, phase, radius);
	153	delta1=10000,delta2=10000;
	154	Double_t d1=1000.,d2=10000.;
b07a036b	155	if (points<=0) return;
51ad6848	156	if (points>0){
	157	phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	158	phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	159	Double_t xd[3],xm[3];
	160	phelix.Evaluate(phase[0][0],xd);
	161	mhelix.Evaluate(phase[0][1],xm);
	162	d1 = (xd[0]-xm[0])(xd[0]-xm[0])+(xd[1]-xm[1])(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]);
	163	}
	164	if (points==2){
	165	phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	166	phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	167	Double_t xd[3],xm[3];
	168	phelix.Evaluate(phase[1][0],xd);
	169	mhelix.Evaluate(phase[1][1],xm);
	170	d2 = (xd[0]-xm[0])(xd[0]-xm[0])+(xd[1]-xm[1])(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]);
	171	}
	172	//
	173	distance2 = TMath::Min(delta1,delta2);
	174	if (delta1<delta2){
	175	//get V0 info
	176	Double_t xd[3],xm[3];
	177	phelix.Evaluate(phase[0][0],xd);
	178	mhelix.Evaluate(phase[0][1], xm);
	179	fXr[0] = 0.5*(xd[0]+xm[0]);
	180	fXr[1] = 0.5*(xd[1]+xm[1]);
	181	fXr[2] = 0.5*(xd[2]+xm[2]);
	182	//
	183	phelix.GetMomentum(phase[0][0],fPP);
	184	mhelix.GetMomentum(phase[0][1],fPM);
	185	phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
	186	fRr = TMath::Sqrt(fXr[0]fXr[0]+fXr[1]fXr[1]);
	187	}
	188	else{
	189	Double_t xd[3],xm[3];
	190	phelix.Evaluate(phase[1][0],xd);
	191	mhelix.Evaluate(phase[1][1], xm);
	192	fXr[0] = 0.5*(xd[0]+xm[0]);
	193	fXr[1] = 0.5*(xd[1]+xm[1]);
	194	fXr[2] = 0.5*(xd[2]+xm[2]);
	195	//
	196	phelix.GetMomentum(phase[1][0], fPP);
	197	mhelix.GetMomentum(phase[1][1], fPM);
	198	phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
	199	fRr = TMath::Sqrt(fXr[0]fXr[0]+fXr[1]fXr[1]);
	200	}
	201	fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
	202	fDist2 = TMath::Sqrt(distance2);
	203	//
	204	//
	205	Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
	206	Float_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
	207	Float_t vnorm2 = v[0]v[0]+v[1]v[1];
	208	Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
	209	vnorm2 = TMath::Sqrt(vnorm2);
	210	Float_t pnorm2 = p[0]p[0]+p[1]p[1];
	211	Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
	212	pnorm2 = TMath::Sqrt(pnorm2);
	213	fPointAngleFi = (v[0]p[0]+v[1]p[1])/(vnorm2*pnorm2);
	214	fPointAngleTh = (v[2]p[2]+vnorm2pnorm2)/(vnorm3*pnorm3);
	215	fPointAngle = (v[0]p[0]+v[1]p[1]+v[2]p[2])/(vnorm3pnorm3);
	216	//
	217	}
	218