]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliGenV0Info.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Several Changes:
[u/mrichter/AliRoot.git] / PWG1 / AliGenV0Info.cxx
CommitLineData
5c7ef659 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///////////////////////////////////////////////////////////////////////////
18/*
19
20Origin: marian.ivanov@cern.ch
21Container classes with MC infomation for V0
22
23
24*/
25
26#if !defined(__CINT__) || defined(__MAKECINT__)
27#include <stdio.h>
28#include <string.h>
29//ROOT includes
30#include "TROOT.h"
31#include "Rtypes.h"
32#include "TFile.h"
33#include "TTree.h"
34#include "TChain.h"
35#include "TCut.h"
36#include "TString.h"
37#include "TStopwatch.h"
38#include "TParticle.h"
39#include "TSystem.h"
40#include "TCanvas.h"
5c7ef659 41#include "TPolyLine3D.h"
42
43//ALIROOT includes
44#include "AliRun.h"
45#include "AliStack.h"
46#include "AliSimDigits.h"
47#include "AliTPCParam.h"
48#include "AliTPC.h"
49#include "AliTPCLoader.h"
50#include "AliDetector.h"
51#include "AliTrackReference.h"
52#include "AliTPCParamSR.h"
53#include "AliTracker.h"
54#include "AliMagF.h"
55#include "AliHelix.h"
56#include "AliTrackPointArray.h"
57
58#endif
59#include "AliGenV0Info.h"
60//
61//
62
63ClassImp(AliGenV0Info)
64
65
66
67
68
69/////////////////////////////////////////////////////////////////////////////////
70AliGenV0Info::AliGenV0Info():
71 fMCd(), //info about daughter particle -
72 fMCm(), //info about mother particle - first particle for V0
73 fMotherP(), //particle info about mother particle
74 fMCDist1(0), //info about closest distance according closest MC - linear DCA
75 fMCDist2(0), //info about closest distance parabolic DCA
76 fMCRr(0), // rec position of the vertex
77 fMCR(0), //exact r position of the vertex
78 fInvMass(0), //reconstructed invariant mass -
79 fPointAngleFi(0), //point angle fi
80 fPointAngleTh(0), //point angle theta
81 fPointAngle(0) //point angle full
82{
83 for (Int_t i=0;i<3; i++){
84 fMCPdr[i]=0;
85 fMCX[i]=0;
86 fMCXr[i]=0;
87 fMCPm[i]=0;
88 fMCAngle[i]=0;
89 fMCPd[i]=0;
90 }
91 fMCPd[3]=0;
92 for (Int_t i=0; i<2;i++){
93 fPdg[i]=0;
94 fLab[i]=0;
95 }
96}
97
98void AliGenV0Info::Update(Float_t vertex[3])
99{
100 //
101 // Update information - calculates derived variables
102 //
103
104 fMCPd[0] = fMCd.GetParticle().Px();
105 fMCPd[1] = fMCd.GetParticle().Py();
106 fMCPd[2] = fMCd.GetParticle().Pz();
107 fMCPd[3] = fMCd.GetParticle().P();
108 //
109 fMCX[0] = fMCd.GetParticle().Vx();
110 fMCX[1] = fMCd.GetParticle().Vy();
111 fMCX[2] = fMCd.GetParticle().Vz();
112 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
113 //
114 fPdg[0] = fMCd.GetParticle().GetPdgCode();
115 fPdg[1] = fMCm.GetParticle().GetPdgCode();
116 //
117 fLab[0] = fMCd.GetParticle().GetUniqueID();
118 fLab[1] = fMCm.GetParticle().GetUniqueID();
119 //
120 //
121 //
122 Double_t x1[3],p1[3];
d390cc7e 123 TParticle& part0 = fMCd.GetParticle();
124 x1[0] = part0.Vx();
125 x1[1] = part0.Vy();
126 x1[2] = part0.Vz();
127 p1[0] = part0.Px();
128 p1[1] = part0.Py();
129 p1[2] = part0.Pz();
130 if (part0.GetPDG()==0) return;
131
132 Double_t sign = (part0.GetPDG()->Charge()>0)? -1:1;
5c7ef659 133 AliHelix dhelix1(x1,p1,sign);
134 //
135 //
136 Double_t x2[3],p2[3];
137 //
d390cc7e 138 TParticle& part1 = fMCm.GetParticle();
139 if (part1.GetPDG()==0) return;
140 x2[0] = part1.Vx();
141 x2[1] = part1.Vy();
142 x2[2] = part1.Vz();
143 p2[0] = part1.Px();
144 p2[1] = part1.Py();
145 p2[2] = part1.Pz();
5c7ef659 146 //
147 //
d390cc7e 148 if (part1.GetPDG()==0) return;
149 sign = (part1.GetPDG()->Charge()>0) ? -1:1;
5c7ef659 150 AliHelix mhelix(x2,p2,sign);
151
152 //
153 //
154 //
155 //find intersection linear
156 //
157 Double_t distance1, distance2;
158 Double_t phase[2][2],radius[2];
159 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
160 Double_t delta1=10000,delta2=10000;
161 if (points>0){
162 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
163 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
164 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
165 }
166 else{
167 fInvMass=-1;
168 return;
169 }
170 if (points==2){
171 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
172 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
173 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
174 }
175 distance1 = TMath::Min(delta1,delta2);
176 //
177 //find intersection parabolic
178 //
179 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
180 delta1=10000,delta2=10000;
181
182 if (points>0){
183 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
184 }
185 if (points==2){
186 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
187 }
188
189 distance2 = TMath::Min(delta1,delta2);
190 //
191 if (delta1<delta2){
192 //get V0 info
193 dhelix1.Evaluate(phase[0][0],fMCXr);
194 dhelix1.GetMomentum(phase[0][0],fMCPdr);
195 mhelix.GetMomentum(phase[0][1],fMCPm);
196 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
197 fMCRr = TMath::Sqrt(radius[0]);
198 }
199 else{
200 dhelix1.Evaluate(phase[1][0],fMCXr);
201 dhelix1.GetMomentum(phase[1][0],fMCPdr);
202 mhelix.GetMomentum(phase[1][1],fMCPm);
203 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
204 fMCRr = TMath::Sqrt(radius[1]);
205 }
206 //
207 //
208 fMCDist1 = TMath::Sqrt(distance1);
209 fMCDist2 = TMath::Sqrt(distance2);
210 //
211 //
212 //
213 Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
214 //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
215 Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
216 Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
217 Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
218 vnorm2 = TMath::Sqrt(vnorm2);
219 Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
220 Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
221 pnorm2 = TMath::Sqrt(pnorm2);
222 //
4c5446dc 223 if (vnorm2>0){
224 fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
225 fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
226 fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
227 }else{
228 fPointAngleFi = 1;
229 fPointAngleTh = 1;
230 fPointAngle = 1;
231 }
5c7ef659 232 Double_t mass1 = fMCd.GetMass();
233 Double_t mass2 = fMCm.GetMass();
234
235
236 Double_t e1 = TMath::Sqrt(mass1*mass1+
237 fMCPd[0]*fMCPd[0]+
238 fMCPd[1]*fMCPd[1]+
239 fMCPd[2]*fMCPd[2]);
240 Double_t e2 = TMath::Sqrt(mass2*mass2+
241 fMCPm[0]*fMCPm[0]+
242 fMCPm[1]*fMCPm[1]+
243 fMCPm[2]*fMCPm[2]);
244
245 fInvMass =
246 (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
247 (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
248 (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
249
250 // fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
251 fInvMass = (e1+e2)*(e1+e2)-fInvMass;
252 if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
253}
254
255
256
257
258
259
260
261