CbmRoot
PairAnalysisPairKF.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 **************************************************************************/
4 
6 // //
7 // PairAnalysis Pair class. Internally it makes use of KFParticle. //
8 // //
10 
11 #include <TDatabasePDG.h>
12 
13 
14 #include "CbmL1.h"
15 #include "CbmL1PFFitter.h"
16 #include "L1Algo.h"
17 #include "L1Field.h"
18 
19 #include "CbmKFParticleInterface.h"
20 #include "CbmKFTrack.h"
21 #include "KFParticle.h"
22 
23 #include "CbmMCTrack.h"
24 #include "CbmVertex.h"
25 #include "PairAnalysisTrack.h"
26 
27 #include "PairAnalysisPairKF.h"
28 
30 
32  : PairAnalysisPair(), fPair(), fD1(), fD2() {
33  //
34  // Default Constructor
35  //
36 }
37 
38 //______________________________________________
40  : PairAnalysisPair(pair), fPair(), fD1(), fD2() {
41  //
42  // Copy Constructor
43  //
45  pair.GetFirstDaughterPid(),
46  pair.GetSecondDaughter(),
47  pair.GetSecondDaughterPid());
48 }
49 
50 //______________________________________________
52  Int_t pid1,
53  PairAnalysisTrack* const particle2,
54  Int_t pid2,
55  Char_t type)
56  : PairAnalysisPair(type), fPair(), fD1(), fD2() {
57  //
58  // Constructor with tracks
59  //
60  SetTracks(particle1, pid1, particle2, pid2);
61 }
62 
63 //______________________________________________
65  //
66  // Default Destructor
67  //
68 }
69 
70 //______________________________________________
72  Int_t pid1,
73  PairAnalysisTrack* const particle2,
74  Int_t pid2) {
75  //
76  // set KF daughters and pair
77  // refParticle1 and 2 are the original tracks. In the case of track rotation
78  // they are needed in the framework
79  //
80  // TODO: think about moving the pid assignement to PairAnalysisTrack and use it here
81  // BUT think about mixed events or LS-pairs
82  // const Double_t mpid1 = TDatabasePDG::Instance()->GetParticle(pid1)->Mass(); (FU) unused
83  // const Double_t mpid2 = TDatabasePDG::Instance()->GetParticle(pid2)->Mass(); (FU) unused
84  const Double_t cpid1 =
85  TDatabasePDG::Instance()->GetParticle(pid1)->Charge() * 3;
86  const Double_t cpid2 =
87  TDatabasePDG::Instance()->GetParticle(pid2)->Charge() * 3;
88 
89  // match charge of track to pid and set mass accordingly
90  fPid1 = pid1;
91  fPid2 = pid2;
92  // Double_t m1 = mpid1; (FU) unused
93  // Double_t m2 = mpid2; (FU) unused
94  if (particle1->Charge() == cpid2) {
95  // m1=mpid2*; (FU) unused
96  fPid1 = pid2;
97  } //TODO: what about 2e-charges
98  if (particle2->Charge() == cpid1) {
99  // m2=mpid1; (FU) unused
100  fPid2 = pid1;
101  }
102 
108  particle1->GetStsTrack(), &fD1, fPid1, kTRUE);
110  particle2->GetStsTrack(), &fD2, fPid2, kTRUE);
111 
112  // references
113  fRefD1 = particle1;
114  fRefD2 = particle2;
115 
116  // build pair
117  fPair.Initialize();
118 
119  fPair.AddDaughter(fD1);
120  fPair.AddDaughter(fD2);
121 
123  // Double_t mass = TDatabasePDG::Instance()->GetParticle(443)->Mass();
124  // Double_t wdth = TDatabasePDG::Instance()->GetParticle(443)->Width();
125  // if(wdth<1.e-6) wdth=mass*0.01; // width<1keV, then 1% mass resolution
126  // fPair.SetMassConstraint( mass, wdth ); //TODO: take from mother pdg code provided to pairanalysis
127 
128  fCharge = (particle1->Charge() * particle2->Charge());
129  fWeight = TMath::Sqrt(particle1->GetWeight() * particle2->GetWeight());
130  // printf("fill pair weight: %.1f * %.1f = %.1f \n",particle1->GetWeight(),particle2->GetWeight(),fWeight);
131 }
132 
133 //______________________________________________
134 void PairAnalysisPairKF::SetMCTracks(const CbmMCTrack* const particle1,
135  const CbmMCTrack* const particle2) {
136  //
137  // build MC pair from daughters
138  // no references are set
139  //
140 
141  //Initialise covariance matrix and set current parameters to 0.0
142  KFParticle kf1;
143  kf1.Initialize();
144  Float_t* par1 = kf1.Parameters();
145  par1[0] = particle1->GetStartX();
146  par1[1] = particle1->GetStartY();
147  par1[2] = particle1->GetStartZ();
148  par1[3] = particle1->GetPx();
149  par1[4] = particle1->GetPy();
150  par1[5] = particle1->GetPz();
151  par1[6] = particle1->GetEnergy();
152  kf1.SetPDG(particle1->GetPdgCode());
153 
154  KFParticle kf2;
155  kf2.Initialize();
156  Float_t* par2 = kf2.Parameters();
157  par2[0] = particle2->GetStartX();
158  par2[1] = particle2->GetStartY();
159  par2[2] = particle2->GetStartZ();
160  par2[3] = particle2->GetPx();
161  par2[4] = particle2->GetPy();
162  par2[5] = particle2->GetPz();
163  par2[6] = particle2->GetEnergy();
164  kf2.SetPDG(particle2->GetPdgCode());
165 
166  // build pair
167  fPair.Initialize();
168  fPair.AddDaughter(kf1);
169  fPair.AddDaughter(kf2);
170 }
171 
172 //______________________________________________
173 void PairAnalysisPairKF::GetThetaPhiCM(Double_t& thetaHE,
174  Double_t& phiHE,
175  Double_t& thetaCS,
176  Double_t& phiCS) const {
177  //
178  // Calculate theta and phi in helicity and Collins-Soper coordinate frame
179  //
180  Double_t px1 = fD1.GetPx();
181  Double_t py1 = fD1.GetPy();
182  Double_t pz1 = fD1.GetPz();
183  Double_t px2 = fD2.GetPx();
184  Double_t py2 = fD2.GetPy();
185  Double_t pz2 = fD2.GetPz();
186  const Double_t d1Mass = TDatabasePDG::Instance()->GetParticle(fPid1)->Mass();
187  const Double_t d2Mass = TDatabasePDG::Instance()->GetParticle(fPid2)->Mass();
188 
189  // first & second daughter 4-mom
190  TLorentzVector p1Mom(
191  px1,
192  py1,
193  pz1,
194  TMath::Sqrt(px1 * px1 + py1 * py1 + pz1 * pz1 + d1Mass * d1Mass));
195  TLorentzVector p2Mom(
196  px2,
197  py2,
198  pz2,
199  TMath::Sqrt(px2 * px2 + py2 * py2 + pz2 * pz2 + d2Mass * d2Mass));
200  // mother 4-momentum vector
201  TLorentzVector motherMom = p1Mom + p2Mom;
202 
204  motherMom, p1Mom, p2Mom, thetaHE, phiHE, thetaCS, phiCS);
205 }
206 
207 
208 //______________________________________________
209 Double_t PairAnalysisPairKF::PsiPair(Double_t /*MagField*/) const {
210  return 0.; /*
211  //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
212  //to ID conversions. Adapted from TRDv0Info class
213  Double_t x, y;//, z;
214  x = fPair.GetX();
215  y = fPair.GetY();
216  // z = fPair.GetZ();
217 
218  Double_t m1[3] = {0,0,0};
219  Double_t m2[3] = {0,0,0};
220 
221  m1[0] = fD1.GetPx();
222  m1[1] = fD1.GetPy();
223  m1[2] = fD1.GetPz();
224 
225  m2[0] = fD2.GetPx();
226  m2[1] = fD2.GetPy();
227  m2[2] = fD2.GetPz();
228 
229  Double_t deltat = 1.;
230  deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
231  TMath::ATan(m1[2]/(TMath::Sqrt(m1[0]*m1[0] + m1[1]*m1[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
232 
233  Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
234 
235  Double_t mom1Prop[3];
236  Double_t mom2Prop[3];
237 
238  ExternalTrackParam *d1 = static_cast<ExternalTrackParam*>(fRefD1.GetObject());
239  ExternalTrackParam *d2 = static_cast<ExternalTrackParam*>(fRefD2.GetObject());
240 
241  ExternalTrackParam nt(*d1), pt(*d2);
242 
243  Double_t fPsiPair = 4.;
244  if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
245  fPsiPair = -5.;
246  if(pt.PropagateTo(radiussum,MagField) == 0)
247  fPsiPair = -5.;
248  pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
249  nt.GetPxPyPz(mom2Prop);
250 
251 
252 
253  Double_t pEle =
254  TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
255  Double_t pPos =
256  TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
257 
258  Double_t scalarproduct =
259  mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
260 
261  Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
262 
263  fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
264 
265  return fPsiPair;
266  */
267 }
268 
269 //______________________________________________
271  //
272  // Calculate the Armenteros-Podolanski Alpha
273  //
274  Int_t qD1 = fD1.GetQ();
275 
276  TVector3 momNeg((qD1 < 0 ? fD1.GetPx() : fD2.GetPx()),
277  (qD1 < 0 ? fD1.GetPy() : fD2.GetPy()),
278  (qD1 < 0 ? fD1.GetPz() : fD2.GetPz()));
279  TVector3 momPos((qD1 < 0 ? fD2.GetPx() : fD1.GetPx()),
280  (qD1 < 0 ? fD2.GetPy() : fD1.GetPy()),
281  (qD1 < 0 ? fD2.GetPz() : fD1.GetPz()));
282  TVector3 momTot(Px(), Py(), Pz());
283 
284  Double_t lQlNeg = momNeg.Dot(momTot) / momTot.Mag();
285  Double_t lQlPos = momPos.Dot(momTot) / momTot.Mag();
286 
287  return ((lQlPos - lQlNeg) / (lQlPos + lQlNeg));
288 }
289 
290 //______________________________________________
292  //
293  // Calculate the Armenteros-Podolanski Pt
294  //
295  Int_t qD1 = fD1.GetQ();
296 
297  TVector3 momNeg((qD1 < 0 ? fD1.GetPx() : fD2.GetPx()),
298  (qD1 < 0 ? fD1.GetPy() : fD2.GetPy()),
299  (qD1 < 0 ? fD1.GetPz() : fD2.GetPz()));
300  TVector3 momTot(Px(), Py(), Pz());
301 
302  return (momNeg.Perp(momTot));
303 }
304 
305 //______________________________________________
306 Double_t PairAnalysisPairKF::PhivPair(Double_t MagField) const {
307  //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
308  //to ID conversions. Angle between ee plane and magnetic field is calculated.
309 
310  //Define local buffer variables for leg properties
311  Double_t px1 = -9999., py1 = -9999., pz1 = -9999.;
312  Double_t px2 = -9999., py2 = -9999., pz2 = -9999.;
313 
314  if (MagField > 0) {
315  if (fD1.GetQ() > 0) {
316  px1 = fD1.GetPx();
317  py1 = fD1.GetPy();
318  pz1 = fD1.GetPz();
319 
320  px2 = fD2.GetPx();
321  py2 = fD2.GetPy();
322  pz2 = fD2.GetPz();
323  } else {
324  px1 = fD2.GetPx();
325  py1 = fD2.GetPy();
326  pz1 = fD2.GetPz();
327 
328  px2 = fD1.GetPx();
329  py2 = fD1.GetPy();
330  pz2 = fD1.GetPz();
331  }
332  } else {
333  if (fD1.GetQ() > 0) {
334  px1 = fD2.GetPx();
335  py1 = fD2.GetPy();
336  pz1 = fD2.GetPz();
337 
338  px2 = fD1.GetPx();
339  py2 = fD1.GetPy();
340  pz2 = fD1.GetPz();
341  } else {
342  px1 = fD1.GetPx();
343  py1 = fD1.GetPy();
344  pz1 = fD1.GetPz();
345 
346  px2 = fD2.GetPx();
347  py2 = fD2.GetPy();
348  pz2 = fD2.GetPz();
349  }
350  }
351 
352  Double_t px = px1 + px2;
353  Double_t py = py1 + py2;
354  Double_t pz = pz1 + pz2;
355  Double_t dppair = TMath::Sqrt(px * px + py * py + pz * pz);
356 
357  //unit vector of (pep+pem)
358  Double_t pl = dppair;
359  Double_t ux = px / pl;
360  Double_t uy = py / pl;
361  Double_t uz = pz / pl;
362  Double_t ax = uy / TMath::Sqrt(ux * ux + uy * uy);
363  Double_t ay = -ux / TMath::Sqrt(ux * ux + uy * uy);
364 
365  //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
366  //definition.
367  //Double_t ptep = iep->Px()*ax + iep->Py()*ay;
368  //Double_t ptem = iem->Px()*ax + iem->Py()*ay;
369 
370  Double_t pxep = px1;
371  Double_t pyep = py1;
372  Double_t pzep = pz1;
373  Double_t pxem = px2;
374  Double_t pyem = py2;
375  Double_t pzem = pz2;
376 
377  //vector product of pep X pem
378  Double_t vpx = pyep * pzem - pzep * pyem;
379  Double_t vpy = pzep * pxem - pxep * pzem;
380  Double_t vpz = pxep * pyem - pyep * pxem;
381  Double_t vp = sqrt(vpx * vpx + vpy * vpy + vpz * vpz);
382  //Double_t thev = acos(vpz/vp);
383 
384  //unit vector of pep X pem
385  Double_t vx = vpx / vp;
386  Double_t vy = vpy / vp;
387  Double_t vz = vpz / vp;
388 
389  //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
390  Double_t wx = uy * vz - uz * vy;
391  Double_t wy = uz * vx - ux * vz;
392  //Double_t wz = ux*vy - uy*vx;
393  //Double_t wl = sqrt(wx*wx+wy*wy+wz*wz);
394  // by construction, (wx,wy,wz) must be a unit vector.
395  // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
396  // should be small if the pair is conversion
397  //
398  Double_t cosPhiV = wx * ax + wy * ay;
399  Double_t phiv = TMath::ACos(cosPhiV);
400 
401  return phiv;
402 }
PairAnalysisPairKF::GetArmAlpha
Double_t GetArmAlpha() const
Definition: PairAnalysisPairKF.cxx:270
PairAnalysisPair::GetFirstDaughterPid
Int_t GetFirstDaughterPid() const
Definition: PairAnalysisPair.h:146
CbmMCTrack::GetStartX
Double_t GetStartX() const
Definition: CbmMCTrack.h:75
PairAnalysisPair::GetFirstDaughter
PairAnalysisTrack * GetFirstDaughter() const
Definition: PairAnalysisPair.h:140
CbmVertex.h
L1Algo.h
PairAnalysisTrack::GetWeight
Double_t GetWeight() const
Definition: PairAnalysisTrack.h:110
CbmKFParticleInterface::SetKFParticleFromStsTrack
static void SetKFParticleFromStsTrack(CbmStsTrack *track, KFParticle *particle, Int_t pdg=211, Bool_t firstPoint=kTRUE)
Definition: CbmKFParticleInterface.cxx:28
L1Field.h
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
PairAnalysisPair::fWeight
Double_t fWeight
Definition: PairAnalysisPair.h:160
PairAnalysisTrack::Charge
Short_t Charge() const
Definition: PairAnalysisTrack.h:107
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
PairAnalysisPairKF::PhivPair
Double_t PhivPair(Double_t MagField) const
Definition: PairAnalysisPairKF.cxx:306
PairAnalysisPair::GetThetaPhiCM
virtual void GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const =0
PairAnalysisPairKF::fPair
KFParticle fPair
Definition: PairAnalysisPairKF.h:143
PairAnalysisPair::fRefD1
TRef fRefD1
beam energy
Definition: PairAnalysisPair.h:164
CbmL1.h
PairAnalysisTrack
Definition: PairAnalysisTrack.h:37
CbmMCTrack::GetPx
Double_t GetPx() const
Definition: CbmMCTrack.h:72
CbmMCTrack::GetPy
Double_t GetPy() const
Definition: CbmMCTrack.h:73
PairAnalysisTrack::GetStsTrack
CbmStsTrack * GetStsTrack() const
Definition: PairAnalysisTrack.h:71
PairAnalysisPair
Definition: PairAnalysisPair.h:25
PairAnalysisPairKF::PairAnalysisPairKF
PairAnalysisPairKF()
CbmMCTrack::GetStartZ
Double_t GetStartZ() const
Definition: CbmMCTrack.h:77
PairAnalysisPairKF::Px
virtual Double_t Px() const
Definition: PairAnalysisPairKF.h:49
PairAnalysisPair::fPid2
Int_t fPid2
Definition: PairAnalysisPair.h:167
PairAnalysisPairKF::Py
virtual Double_t Py() const
Definition: PairAnalysisPairKF.h:50
PairAnalysisPairKF::GetArmPt
Double_t GetArmPt() const
Definition: PairAnalysisPairKF.cxx:291
PairAnalysisPair::fCharge
Short_t fCharge
Definition: PairAnalysisPair.h:158
PairAnalysisPairKF::GetThetaPhiCM
void GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
Definition: PairAnalysisPairKF.cxx:173
CbmL1PFFitter.h
CbmKFTrack.h
CbmKFParticleInterface.h
PairAnalysisPairKF::fD2
KFParticle fD2
Definition: PairAnalysisPairKF.h:145
ClassImp
ClassImp(PairAnalysisPairKF) PairAnalysisPairKF
Definition: PairAnalysisPairKF.cxx:29
PairAnalysisTrack.h
CbmMCTrack::GetStartY
Double_t GetStartY() const
Definition: CbmMCTrack.h:76
PairAnalysisPair::fPid1
Int_t fPid1
Definition: PairAnalysisPair.h:166
PairAnalysisPairKF::SetMCTracks
void SetMCTracks(const CbmMCTrack *const particle1, const CbmMCTrack *const particle2)
Definition: PairAnalysisPairKF.cxx:134
CbmMCTrack.h
PairAnalysisPair::GetSecondDaughterPid
Int_t GetSecondDaughterPid() const
Definition: PairAnalysisPair.h:147
PairAnalysisPairKF::Pz
virtual Double_t Pz() const
Definition: PairAnalysisPairKF.h:51
PairAnalysisPair::fRefD2
TRef fRefD2
Definition: PairAnalysisPair.h:165
PairAnalysisPairKF::fD1
KFParticle fD1
Definition: PairAnalysisPairKF.h:144
CbmMCTrack
Definition: CbmMCTrack.h:34
PairAnalysisPairKF::~PairAnalysisPairKF
virtual ~PairAnalysisPairKF()
Definition: PairAnalysisPairKF.cxx:64
PairAnalysisPairKF::SetTracks
void SetTracks(PairAnalysisTrack *const particle1, Int_t pid1, PairAnalysisTrack *const particle2, Int_t pid2)
Definition: PairAnalysisPairKF.cxx:71
CbmMCTrack::GetEnergy
Double_t GetEnergy() const
Definition: CbmMCTrack.h:165
PairAnalysisPairKF
Definition: PairAnalysisPairKF.h:28
PairAnalysisPairKF::PsiPair
Double_t PsiPair(Double_t MagField) const
Definition: PairAnalysisPairKF.cxx:209
CbmMCTrack::GetPz
Double_t GetPz() const
Definition: CbmMCTrack.h:74
PairAnalysisPair::GetSecondDaughter
PairAnalysisTrack * GetSecondDaughter() const
Definition: PairAnalysisPair.h:143
PairAnalysisPairKF.h