[u/mrichter/AliRoot.git] / STEER / AliCascadeVertexer.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.                  *
 **************************************************************************/

//-------------------------------------------------------------------------
//               Implementation of the cascade vertexer class
//          Reads V0s and tracks, writes out cascade vertices
//                     Fills the ESD with the cascades 
//    Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
//-------------------------------------------------------------------------
#include <TObjArray.h>
#include <TTree.h>

#include "AliESD.h"
#include "AliESDv0.h"
#include "AliESDcascade.h"
#include "AliCascadeVertexer.h"

ClassImp(AliCascadeVertexer)

Int_t AliCascadeVertexer::V0sTracks2CascadeVertices(AliESD *event) {
  //--------------------------------------------------------------------
  // This function reconstructs cascade vertices
  //      Adapted to the ESD by I.Belikov (Jouri.Belikov@cern.ch)
  //--------------------------------------------------------------------
   Double_t b=event->GetMagneticField();
   Int_t nV0=(Int_t)event->GetNumberOfV0s();
   TObjArray vtcs(nV0);
   Int_t i;
   for (i=0; i<nV0; i++) {
       AliESDv0 *v=event->GetV0(i);
       if (v->GetD(fX,fY,fZ)<fDV0min) continue;
       vtcs.AddLast(v);
   }
   nV0=vtcs.GetEntriesFast();

   Int_t nentr=(Int_t)event->GetNumberOfTracks();
   TArrayI trk(nentr); Int_t ntr=0;
   for (i=0; i<nentr; i++) {
       AliESDtrack *esdtr=event->GetTrack(i);
       UInt_t status=esdtr->GetStatus();
       UInt_t flags=AliESDtrack::kITSin|AliESDtrack::kTPCin|
                    AliESDtrack::kTPCpid|AliESDtrack::kESDpid;

       if ((status&AliESDtrack::kITSrefit)==0)
          if (flags!=status) continue;

       if (TMath::Abs(esdtr->GetD(fX,fY,b))<fDBachMin) continue;

       trk[ntr++]=i;
   }   

   Double_t massLambda=1.11568;
   Int_t ncasc=0;

   // Looking for the cascades...
   for (i=0; i<nV0; i++) {
      AliESDv0 *v=(AliESDv0*)vtcs.UncheckedAt(i);
      v->ChangeMassHypothesis(kLambda0); // the v0 must be Lambda 
      if (TMath::Abs(v->GetEffMass()-massLambda)>fMassWin) continue; 
      for (Int_t j=0; j<ntr; j++) {
	 Int_t bidx=trk[j];
	 AliESDtrack *btrk=event->GetTrack(bidx);

         if (btrk->GetSign()>0) continue;  // bachelor's charge 
          
	 AliESDv0 v0(*v), *pv0=&v0;
         AliExternalTrackParam bt(*btrk), *pbt=&bt;

         Double_t dca=PropagateToDCA(pv0,pbt,b);
         if (dca > fDCAmax) continue;

         AliESDcascade cascade(*pv0,*pbt,bidx);
         if (cascade.GetChi2() > fChi2max) continue;

	 Double_t x,y,z; cascade.GetXYZ(x,y,z); 
         Double_t r2=x*x + y*y; 
         if (r2 > fRmax*fRmax) continue;   // condition on fiducial zone
         if (r2 < fRmin*fRmin) continue;

         {
         Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
         if (r2 > (x1*x1+y1*y1)) continue;
         //if ((z-fZ)*(z-fZ) > (z1-fZ)*(z1-fZ)) continue;
         }

	 Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz);
         Double_t p2=px*px+py*py+pz*pz;
         Double_t cost=((x-fX)*px + (y-fY)*py + (z-fZ)*pz)/
               TMath::Sqrt(p2*((x-fX)*(x-fX) + (y-fY)*(y-fY) + (z-fZ)*(z-fZ)));
        if (cost<fCPAmax) continue; //condition on the cascade pointing angle 
         //cascade.ChangeMassHypothesis(); //default is Xi

         event->AddCascade(&cascade);

         ncasc++;

      }
   }

   // Looking for the anti-cascades...
   for (i=0; i<nV0; i++) {
      AliESDv0 *v=(AliESDv0*)vtcs.UncheckedAt(i);
      v->ChangeMassHypothesis(kLambda0Bar); //the v0 must be anti-Lambda 
      if (TMath::Abs(v->GetEffMass()-massLambda)>fMassWin) continue; 
      for (Int_t j=0; j<ntr; j++) {
	 Int_t bidx=trk[j];
	 AliESDtrack *btrk=event->GetTrack(bidx);

         if (btrk->GetSign()<0) continue;  // bachelor's charge 
          
	 AliESDv0 v0(*v), *pv0=&v0;
         AliESDtrack bt(*btrk), *pbt=&bt;

         Double_t dca=PropagateToDCA(pv0,pbt,b);
         if (dca > fDCAmax) continue;

         AliESDcascade cascade(*pv0,*pbt,bidx);
         if (cascade.GetChi2() > fChi2max) continue;

	 Double_t x,y,z; cascade.GetXYZ(x,y,z); 
         Double_t r2=x*x + y*y; 
         if (r2 > fRmax*fRmax) continue;   // condition on fiducial zone
         if (r2 < fRmin*fRmin) continue;

         {
         Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
         if (r2 > (x1*x1+y1*y1)) continue;
         if (z*z > z1*z1) continue;
         }

	 Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz);
         Double_t p2=px*px+py*py+pz*pz;
         Double_t cost=((x-fX)*px + (y-fY)*py + (z-fZ)*pz)/
               TMath::Sqrt(p2*((x-fX)*(x-fX) + (y-fY)*(y-fY) + (z-fZ)*(z-fZ)));

         if (cost<fCPAmax) continue; //condition on the cascade pointing angle 
         //cascade.ChangeMassHypothesis(); //default is Xi

         event->AddCascade(&cascade);

         ncasc++;

      }
   }

Info("V0sTracks2CascadeVertices","Number of reconstructed cascades: %d",ncasc);

   return 0;
}


inline Double_t det(Double_t a00, Double_t a01, Double_t a10, Double_t a11){
  // determinant 2x2
  return a00*a11 - a01*a10;
}

inline Double_t det (Double_t a00,Double_t a01,Double_t a02,
                     Double_t a10,Double_t a11,Double_t a12,
                     Double_t a20,Double_t a21,Double_t a22) {
  // determinant 3x3
  return 
  a00*det(a11,a12,a21,a22)-a01*det(a10,a12,a20,a22)+a02*det(a10,a11,a20,a21);
}


Double_t AliCascadeVertexer::
PropagateToDCA(AliESDv0 *v, AliExternalTrackParam *t, Double_t b) {
  //--------------------------------------------------------------------
  // This function returns the DCA between the V0 and the track
  //--------------------------------------------------------------------
  Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
  Double_t r[3]; t->GetXYZ(r);
  Double_t x1=r[0], y1=r[1], z1=r[2];
  Double_t p[3]; t->GetPxPyPz(p);
  Double_t px1=p[0], py1=p[1], pz1=p[2];
  
  Double_t x2,y2,z2;     // position and momentum of V0
  Double_t px2,py2,pz2;
  
  v->GetXYZ(x2,y2,z2);
  v->GetPxPyPz(px2,py2,pz2);
 
// calculation dca
   
  Double_t dd= det(x2-x1,y2-y1,z2-z1,px1,py1,pz1,px2,py2,pz2);
  Double_t ax= det(py1,pz1,py2,pz2);
  Double_t ay=-det(px1,pz1,px2,pz2);
  Double_t az= det(px1,py1,px2,py2);

  Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az);

//points of the DCA
  Double_t t1 = det(x2-x1,y2-y1,z2-z1,px2,py2,pz2,ax,ay,az)/
                det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
  
  x1 += px1*t1; y1 += py1*t1; //z1 += pz1*t1;
  

  //propagate track to the points of DCA

  x1=x1*cs1 + y1*sn1;
  if (!t->PropagateTo(x1,b)) {
    Error("PropagateToDCA","Propagation failed !");
    return 1.e+33;
  }  

  return dca;
}
Commit	Line	Data
c7bafca9	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	//-------------------------------------------------------------------------
	17	// Implementation of the cascade vertexer class
	18	// Reads V0s and tracks, writes out cascade vertices
	19	// Fills the ESD with the cascades
	20	// Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
	21	//-------------------------------------------------------------------------
	22	#include <TObjArray.h>
	23	#include <TTree.h>
	24
	25	#include "AliESD.h"
	26	#include "AliESDv0.h"
	27	#include "AliESDcascade.h"
	28	#include "AliCascadeVertexer.h"
	29
	30	ClassImp(AliCascadeVertexer)
	31
	32	Int_t AliCascadeVertexer::V0sTracks2CascadeVertices(AliESD *event) {
	33	//--------------------------------------------------------------------
	34	// This function reconstructs cascade vertices
	35	// Adapted to the ESD by I.Belikov (Jouri.Belikov@cern.ch)
	36	//--------------------------------------------------------------------
	37	Double_t b=event->GetMagneticField();
	38	Int_t nV0=(Int_t)event->GetNumberOfV0s();
	39	TObjArray vtcs(nV0);
	40	Int_t i;
	41	for (i=0; i<nV0; i++) {
	42	AliESDv0 *v=event->GetV0(i);
	43	if (v->GetD(fX,fY,fZ)<fDV0min) continue;
	44	vtcs.AddLast(v);
	45	}
	46	nV0=vtcs.GetEntriesFast();
	47
	48	Int_t nentr=(Int_t)event->GetNumberOfTracks();
	49	TArrayI trk(nentr); Int_t ntr=0;
	50	for (i=0; i<nentr; i++) {
	51	AliESDtrack *esdtr=event->GetTrack(i);
	52	UInt_t status=esdtr->GetStatus();
	53	UInt_t flags=AliESDtrack::kITSin\|AliESDtrack::kTPCin\|
	54	AliESDtrack::kTPCpid\|AliESDtrack::kESDpid;
	55
	56	if ((status&AliESDtrack::kITSrefit)==0)
	57	if (flags!=status) continue;
	58
	59	if (TMath::Abs(esdtr->GetD(fX,fY,b))<fDBachMin) continue;
	60
	61	trk[ntr++]=i;
	62	}
	63
	64	Double_t massLambda=1.11568;
65	Int_t ncasc=0;
66
67	// Looking for the cascades...
68	for (i=0; i<nV0; i++) {
69	AliESDv0 v=(AliESDv0)vtcs.UncheckedAt(i);
70	v->ChangeMassHypothesis(kLambda0); // the v0 must be Lambda
71	if (TMath::Abs(v->GetEffMass()-massLambda)>fMassWin) continue;
72	for (Int_t j=0; j<ntr; j++) {
73	Int_t bidx=trk[j];
74	AliESDtrack *btrk=event->GetTrack(bidx);
75
76	if (btrk->GetSign()>0) continue; // bachelor's charge
77
78	AliESDv0 v0(v), pv0=&v0;
79	AliExternalTrackParam bt(btrk), pbt=&bt;
80
81	Double_t dca=PropagateToDCA(pv0,pbt,b);
82	if (dca > fDCAmax) continue;
83
84	AliESDcascade cascade(pv0,pbt,bidx);
85	if (cascade.GetChi2() > fChi2max) continue;
86
87	Double_t x,y,z; cascade.GetXYZ(x,y,z);
88	Double_t r2=xx + yy;
89	if (r2 > fRmax*fRmax) continue; // condition on fiducial zone
90	if (r2 < fRmin*fRmin) continue;
91
92	{
93	Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
94	if (r2 > (x1x1+y1y1)) continue;
95	//if ((z-fZ)(z-fZ) > (z1-fZ)(z1-fZ)) continue;
96	}
97
98	Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz);
99	Double_t p2=pxpx+pypy+pz*pz;
100	Double_t cost=((x-fX)px + (y-fY)py + (z-fZ)*pz)/
101	TMath::Sqrt(p2((x-fX)(x-fX) + (y-fY)(y-fY) + (z-fZ)(z-fZ)));
102	if (cost<fCPAmax) continue; //condition on the cascade pointing angle
103	//cascade.ChangeMassHypothesis(); //default is Xi
104
105	event->AddCascade(&cascade);
106
107	ncasc++;
108
109	}
110	}
111
112	// Looking for the anti-cascades...
113	for (i=0; i<nV0; i++) {
114	AliESDv0 v=(AliESDv0)vtcs.UncheckedAt(i);
115	v->ChangeMassHypothesis(kLambda0Bar); //the v0 must be anti-Lambda
116	if (TMath::Abs(v->GetEffMass()-massLambda)>fMassWin) continue;
117	for (Int_t j=0; j<ntr; j++) {
118	Int_t bidx=trk[j];
119	AliESDtrack *btrk=event->GetTrack(bidx);
120
121	if (btrk->GetSign()<0) continue; // bachelor's charge
122
123	AliESDv0 v0(v), pv0=&v0;
124	AliESDtrack bt(btrk), pbt=&bt;
125
126	Double_t dca=PropagateToDCA(pv0,pbt,b);
127	if (dca > fDCAmax) continue;
128
129	AliESDcascade cascade(pv0,pbt,bidx);
130	if (cascade.GetChi2() > fChi2max) continue;
131
132	Double_t x,y,z; cascade.GetXYZ(x,y,z);
133	Double_t r2=xx + yy;
134	if (r2 > fRmax*fRmax) continue; // condition on fiducial zone
135	if (r2 < fRmin*fRmin) continue;
136
137	{
138	Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
139	if (r2 > (x1x1+y1y1)) continue;
140	if (zz > z1z1) continue;
141	}
142
143	Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz);
144	Double_t p2=pxpx+pypy+pz*pz;
145	Double_t cost=((x-fX)px + (y-fY)py + (z-fZ)*pz)/
146	TMath::Sqrt(p2((x-fX)(x-fX) + (y-fY)(y-fY) + (z-fZ)(z-fZ)));
147
148	if (cost<fCPAmax) continue; //condition on the cascade pointing angle
149	//cascade.ChangeMassHypothesis(); //default is Xi
150
151	event->AddCascade(&cascade);
152
153	ncasc++;
154
155	}
156	}
157
158	Info("V0sTracks2CascadeVertices","Number of reconstructed cascades: %d",ncasc);
159
160	return 0;
161	}
162
163
164	inline Double_t det(Double_t a00, Double_t a01, Double_t a10, Double_t a11){
165	// determinant 2x2
166	return a00a11 - a01a10;
167	}
168
169	inline Double_t det (Double_t a00,Double_t a01,Double_t a02,
170	Double_t a10,Double_t a11,Double_t a12,
171	Double_t a20,Double_t a21,Double_t a22) {
172	// determinant 3x3
173	return
174	a00det(a11,a12,a21,a22)-a01det(a10,a12,a20,a22)+a02*det(a10,a11,a20,a21);
175	}
176
177
178
179
180	Double_t AliCascadeVertexer::
181	PropagateToDCA(AliESDv0 v, AliExternalTrackParam t, Double_t b) {
182	//--------------------------------------------------------------------
183	// This function returns the DCA between the V0 and the track
184	//--------------------------------------------------------------------
185	Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
186	Double_t r[3]; t->GetXYZ(r);
187	Double_t x1=r[0], y1=r[1], z1=r[2];
188	Double_t p[3]; t->GetPxPyPz(p);
189	Double_t px1=p[0], py1=p[1], pz1=p[2];
190
191	Double_t x2,y2,z2; // position and momentum of V0
192	Double_t px2,py2,pz2;
193
194	v->GetXYZ(x2,y2,z2);
195	v->GetPxPyPz(px2,py2,pz2);
196
197	// calculation dca
198
199	Double_t dd= det(x2-x1,y2-y1,z2-z1,px1,py1,pz1,px2,py2,pz2);
200	Double_t ax= det(py1,pz1,py2,pz2);
201	Double_t ay=-det(px1,pz1,px2,pz2);
202	Double_t az= det(px1,py1,px2,py2);
203
204	Double_t dca=TMath::Abs(dd)/TMath::Sqrt(axax + ayay + az*az);
205
206	//points of the DCA
207	Double_t t1 = det(x2-x1,y2-y1,z2-z1,px2,py2,pz2,ax,ay,az)/
208	det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
209
210	x1 += px1t1; y1 += py1t1; //z1 += pz1*t1;
211
212
213	//propagate track to the points of DCA
214
215	x1=x1cs1 + y1sn1;
216	if (!t->PropagateTo(x1,b)) {
217	Error("PropagateToDCA","Propagation failed !");
218	return 1.e+33;
219	}
220
221	return dca;
222	}
223
224
225
226
227
228
229
230
231
232
233