don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliGenKinkInfo.cxx
CommitLineData
7cc34f08 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
22
23*/
24
25#if !defined(__CINT__) || defined(__MAKECINT__)
26#include <stdio.h>
27#include <string.h>
28//ROOT includes
29#include "TROOT.h"
30#include "Rtypes.h"
31#include "TFile.h"
32#include "TTree.h"
33#include "TChain.h"
34#include "TCut.h"
35#include "TString.h"
36#include "TStopwatch.h"
37#include "TParticle.h"
38#include "TDatabasePDG.h"
39#include "TSystem.h"
40#include "TCanvas.h"
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 "AliGenKinkInfo.h"
60
61//
62//
63
64ClassImp(AliGenKinkInfo)
65
66
67
68/////////////////////////////////////////////////////////////////////////////////
69
70AliGenKinkInfo::AliGenKinkInfo():
71 fMCd(), //info about daughter particle - second particle for V0
72 fMCm(), //info about mother particle - first particle for V0
73 fMCDist1(0), //info about closest distance according closest MC - linear DCA
74 fMCDist2(0), //info about closest distance parabolic DCA
75 fMCRr(0), // rec position of the vertex
76 fMCR(0) //exact r position of the vertex
77{
78 //
79 // default constructor
80 //
81 for (Int_t i=0;i<3;i++){
82 fMCPdr[i]=0;
83 fMCPd[i]=0;
84 fMCX[i]=0;
85 fMCPm[i]=0;
86 fMCAngle[i]=0;
4a9220d4 87 fMCXr[i] = 0;
7cc34f08 88 }
89 for (Int_t i=0; i<2; i++) {
90 fPdg[i]= 0;
91 fLab[i]=0;
92 }
93}
94
95void AliGenKinkInfo::Update()
96{
97 //
98 // Update information
99 // some redundancy - faster acces to the values in analysis code
100 //
101 fMCPd[0] = fMCd.GetParticle().Px();
102 fMCPd[1] = fMCd.GetParticle().Py();
103 fMCPd[2] = fMCd.GetParticle().Pz();
104 fMCPd[3] = fMCd.GetParticle().P();
105 //
106 fMCX[0] = fMCd.GetParticle().Vx();
107 fMCX[1] = fMCd.GetParticle().Vy();
108 fMCX[2] = fMCd.GetParticle().Vz();
109 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
110 //
111 fPdg[0] = fMCd.GetParticle().GetPdgCode();
112 fPdg[1] = fMCm.GetParticle().GetPdgCode();
113 //
114 fLab[0] = fMCd.GetParticle().GetUniqueID();
115 fLab[1] = fMCm.GetParticle().GetUniqueID();
116 //
117 //
118 //
119 Double_t x1[3],p1[3];
120 TParticle& pdaughter = fMCd.GetParticle();
121 x1[0] = pdaughter.Vx();
122 x1[1] = pdaughter.Vy();
123 x1[2] = pdaughter.Vz();
124 p1[0] = pdaughter.Px();
125 p1[1] = pdaughter.Py();
126 p1[2] = pdaughter.Pz();
cb474242 127 Double_t sign = 0;
128 sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
7cc34f08 129 AliHelix dhelix1(x1,p1,sign);
130 //
131 //
132 Double_t x2[3],p2[3];
133 //
134 TParticle& pmother = fMCm.GetParticle();
135 x2[0] = pmother.Vx();
136 x2[1] = pmother.Vy();
137 x2[2] = pmother.Vz();
138 p2[0] = pmother.Px();
139 p2[1] = pmother.Py();
140 p2[2] = pmother.Pz();
141 //
142 const AliTrackReference & pdecay = fMCm.GetTRdecay();
143 x2[0] = pdecay.X();
144 x2[1] = pdecay.Y();
145 x2[2] = pdecay.Z();
146 p2[0] = pdecay.Px();
147 p2[1] = pdecay.Py();
148 p2[2] = pdecay.Pz();
149 //
150 sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
151 AliHelix mhelix(x2,p2,sign);
152
153 //
154 //
155 //
156 //find intersection linear
157 //
158 Double_t distance1, distance2;
4a9220d4 159 Double_t phase[2][2] = { {0,0},{0,0} };
160 Double_t radius[2] = {0};
7cc34f08 161 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
162 Double_t delta1=10000,delta2=10000;
163 if (points>0){
164 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
165 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
166 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
167 }
168 if (points==2){
169 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
170 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
171 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
172 }
173 distance1 = TMath::Min(delta1,delta2);
174 //
175 //find intersection parabolic
176 //
177 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
178 delta1=10000,delta2=10000;
179
180 if (points>0){
181 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
182 }
183 if (points==2){
184 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
185 }
186
187 distance2 = TMath::Min(delta1,delta2);
188 //
189 if (delta1<delta2){
190 //get V0 info
191 dhelix1.Evaluate(phase[0][0],fMCXr);
192 dhelix1.GetMomentum(phase[0][0],fMCPdr);
193 mhelix.GetMomentum(phase[0][1],fMCPm);
194 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
195 fMCRr = TMath::Sqrt(radius[0]);
196 }
197 else{
198 dhelix1.Evaluate(phase[1][0],fMCXr);
199 dhelix1.GetMomentum(phase[1][0],fMCPdr);
200 mhelix.GetMomentum(phase[1][1],fMCPm);
201 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
202 fMCRr = TMath::Sqrt(radius[1]);
203 }
204 //
205 //
206 fMCDist1 = TMath::Sqrt(distance1);
207 fMCDist2 = TMath::Sqrt(distance2);
208
209}
210
211
212Float_t AliGenKinkInfo::GetQt(){
213 //
214 //
215 Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
216 return TMath::Sin(fMCAngle[2])*momentumd;
217}
218
219