StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StXiFinderMaker.cxx
1 //
7 // Cuts can be found in the code by comments beginning with "Cut:"
8 //
9 //
10 
11 
12 #include "StXiFinderMaker.h"
13 #include "StMessMgr.h"
14 #include "StEvent/StEventTypes.h"
15 #include "TMath.h"
16 #include "TVector2.h"
17 #include "tables/St_exi_exipar_Table.h"
18 #include "PhysicalConstants.h"
19 #include "phys_constants.h"
20 
21 
22 StSPtrVecXiVertex* vecXi=0;
23 
24 
25 
26 //_____________________________________________________________________________
27 StXiFinderMaker::StXiFinderMaker(const char *name):StV0FinderMaker(name),
28 exipar(0),parsXi(0),xiVertex(0),det_id_xi(0){}
29 
30 //_____________________________________________________________________________
31 StXiFinderMaker::~StXiFinderMaker() {}
32 //_____________________________________________________________________________
33 Int_t StXiFinderMaker::Init()
34 {bool a,b,c;
35  if ((useTracker!=kTrackerUseTPT) && (useTracker!=kTrackerUseITTF) && (useTracker!=kTrackerUseBOTH))
36  {gMessMgr->Error("StXiFinderMaker::Init() : wrong TrackerUsage parameter set.");
37  return kStErr;
38  }
39  if ((useSVT!=kNoSVT) && (useSVT!=kUseSVT))
40  {gMessMgr->Error("StXiFinderMaker::Init() : wrong SVTUsage parameter set.");
41  return kStErr;
42  }
43  if ((useEventModel!=kUseStEvent) && (useEventModel!=kUseMuDst))
44  {gMessMgr->Error("StXiFinderMaker::Init() : wrong EventModelUsage parameter set.");
45  return kStErr;
46  }
47  if ((useLikesign!=kLikesignUseStandard) && (useLikesign!=kLikesignUseLikesign))
48  {gMessMgr->Error("StXiFinderMaker::Init() : wrong LikesignUsage parameter set.");
49  return kStErr;
50  }
51  if ((useRotating!=kRotatingUseStandard) && (useRotating!=kRotatingUseRotating) && (useRotating!=kRotatingUseSymmetry) && (useRotating!=kRotatingUseRotatingAndSymmetry))
52  {gMessMgr->Error("StXiFinderMaker::Init() : wrong RotatingUsage parameter set.");
53  return kStErr;
54  }
55 
56  if (useTracker == kTrackerUseTPT) gMessMgr->Info()<<"StXiFinderMaker : use TPT tracks."<<endm;
57  if (useTracker == kTrackerUseITTF) gMessMgr->Info()<<"StXiFinderMaker : use ITTF tracks."<<endm;
58  if (useTracker == kTrackerUseBOTH) gMessMgr->Info()<<"StXiFinderMaker : use TPT *and* ITTF tracks."<<endm;
59  if (useSVT == kUseSVT) gMessMgr->Info()<<"StXiFinderMaker : use SVT points if possible."<<endm;
60  if (useSVT == kNoSVT) gMessMgr->Info()<<"StXiFinderMaker : do not use SVT points."<<endm;
61  if (useEventModel == kUseStEvent) gMessMgr->Info()<<"StXiFinderMaker : expect StEvent files in input."<<endm;
62  if (useEventModel == kUseMuDst) gMessMgr->Info()<<"StXiFinderMaker : expect MuDst files in input."<<endm;
63  if (useLikesign == kLikesignUseLikesign) gMessMgr->Info()<<"StXiFinderMaker : does like-sign finding."<<endm;
64  if (useRotating == kRotatingUseRotating) gMessMgr->Info()<<"StXiFinderMaker : does rotating finding."<<endm;
65  if (useRotating == kRotatingUseSymmetry) gMessMgr->Info()<<"StXiFinderMaker : does symmetry finding."<<endm;
66  if (useRotating == kRotatingUseRotatingAndSymmetry) gMessMgr->Info()<<"StXiFinderMaker : does rotating + symmetry finding."<<endm;
67 
68  if (useLanguage != kLanguageUseSpecial)
69  {a=(bool)(1&(useLanguage>>2));
70  b=(bool)(1&(useLanguage>>1));
71  c=(bool)(1&useLanguage);
72  useV0Language=2*(!(a^c))+(a|c);
73  useXiLanguage=4*(b&(!(a^c)))+2*(a&b&(!c))+(a|c);
74  }
75  switch (useLanguage)
76  {case kLanguageUseOldRun : gMessMgr->Info()<<"StXiFinderMaker : Fortran run."<<endm;
77  break;
78  case kLanguageUseRun : gMessMgr->Info()<<"StXiFinderMaker : C++ run."<<endm;
79  break;
80  case kLanguageUseTestV0Finder : gMessMgr->Info()<<"StXiFinderMaker : Test V0Finder."<<endm;
81  break;
82  case kLanguageUseTestXiFinder : gMessMgr->Info()<<"StXiFinderMaker : Test XiFinder."<<endm;
83  break;
84  case kLanguageUseTestBothFinders : gMessMgr->Info()<<"StXiFinderMaker : Test V0Finder and XiFinder."<<endm;
85  break;
86  case kLanguageUseSpecial : break;
87  default : gMessMgr->Error("StXiFinderMaker::Init() : wrong LanguageUsage parameter set.");
88  return kStErr;
89  }
90  if ((useXiLanguage!=kXiLanguageUseFortran) && (useXiLanguage!=kXiLanguageUseCppOnFortranV0) && (useXiLanguage!=kXiLanguageUseCppOnCppV0) && (useXiLanguage!=kXiLanguageUseFortranAndCppOnFortranV0) && (useXiLanguage!=kXiLanguageUseFortranAndCppOnCppV0) && (useXiLanguage!=kXiLanguageUseBothCpp) && (useXiLanguage!=kXiLanguageUseAll))
91  {gMessMgr->Error("StXiFinderMaker::Init() : wrong XiLanguageUsage parameter set.");
92  return kStErr;
93  }
94  switch (useV0Language)
95  {case kV0LanguageUseFortran : if ((useXiLanguage!=kXiLanguageUseFortran) && (useXiLanguage!=kXiLanguageUseCppOnFortranV0) && (useXiLanguage!=kXiLanguageUseFortranAndCppOnFortranV0))
96  {gMessMgr->Info()<<"StXiFinderMaker : BE CAREFUL : impossible combination asked."<<endm;
97  gMessMgr->Info()<<"StXiFinderMaker : Set it to testXiFinder."<<endm;
98  useXiLanguage=kXiLanguageUseFortranAndCppOnFortranV0;
99  }
100  break;
101  case kV0LanguageUseCpp : if (useXiLanguage!=kXiLanguageUseCppOnCppV0)
102  {gMessMgr->Info()<<"StXiFinderMaker : BE CAREFUL : impossible combination asked."<<endm;
103  gMessMgr->Info()<<"StXiFinderMaker : Set it to normalRun."<<endm;
104  useXiLanguage=kXiLanguageUseCppOnCppV0;
105  }
106  break;
107  case kV0LanguageUseBoth : break;
108  default : gMessMgr->Error("StXiFinderMaker::Init() : wrong V0LanguageUsage parameter set.");
109  return kStErr;
110  }
111  if (1&useV0Language) gMessMgr->Info()<<"StXiFinderMaker : Will store Fortran V0s."<<endm;
112  if (2&useV0Language) gMessMgr->Info()<<"StXiFinderMaker : Will store C++ V0s."<<endm;
113  if (1&useXiLanguage) gMessMgr->Info()<<"StXiFinderMaker : Will store Fortran Xis."<<endm;
114  if (2&useXiLanguage) gMessMgr->Info()<<"StXiFinderMaker : Will store C++ Xis made with Fortran V0s."<<endm;
115  if (4&useXiLanguage) gMessMgr->Info()<<"StXiFinderMaker : Will store C++ Xis made with C++ V0s."<<endm;
116 
117  if (useEventModel) //initialize mMuDstMaker
118  {mMuDstMaker = (StMuDstMaker*)GetMaker("myMuDstMaker");
119  if(!mMuDstMaker) gMessMgr->Warning("StXiFinderMaker::Init can't find a valid MuDst.");
120  }
121 
122  return StMaker::Init();
123  }
124 //=============================================================================
125 Int_t StXiFinderMaker::InitRun(int runumber) {
126  exipar = (St_exi_exipar*) GetDataBase("global/vertices/exipar");
127  if (!exipar)
128  {gMessMgr->Error("StXiFinderMaker::Init() : could not find exipar in database.");
129  return kStErr;
130  }
131  return StMaker::InitRun(runumber);
132  }
133 //=============================================================================
134 //_____________________________________________________________________________
135 Int_t StXiFinderMaker::Make() {
136 
137  Int_t iRes;
138 
139  // Prepare event and track variables
140  iRes = Prepare();
141  if (iRes != kStOk) return iRes;
142 
143  StSPtrVecXiVertex& xiVertices = event->xiVertices();
144  gMessMgr->Info()<<"StXiFinderMaker : coming in I have "<<xiVertices.size()<<" Xis."<<endm;
157  if (!(1&useXiLanguage) && !((4&useXiLanguage) && (!(1&useV0Language))))
158  {// Erase existing Xis
159  gMessMgr->Info() << "StXiFinderMaker : pre-existing Xis deleted." << endm;
160  StSPtrVecXiVertex xiVertices2;
161  xiVertices = xiVertices2;
162  }
163  vecXi = &xiVertices;
164 
165  // Call the V0-finding, which will in turn call
166  // the UseV0() member function for each V0 found
167 
168  if (2&useXiLanguage)
169  {StSPtrVecV0Vertex& v0Vertices = event->v0Vertices();
170  unsigned int nV0s = v0Vertices.size();
171  det_id_v0 = 1; // for lack of any further information
172  for (unsigned int i=0; i<nV0s; i++)
173  {v0Vertex = v0Vertices[i];
174  if (v0Vertex) UseV0();
175  }
176  }
177  if ((4&useXiLanguage) || (2&useV0Language))
178  {iRes = StV0FinderMaker::Make();
179  if (iRes != kStOk) return iRes;
180  }
181 
182  gMessMgr->Info()<<"StXiFinderMaker : now I have "<<xiVertices.size()<<" Xis."<<endm;
183 
184  return kStOk;
185 }
186 //_____________________________________________________________________________
187 Bool_t StXiFinderMaker::UseV0() {
188 
189  Bool_t usedV0 = kFALSE;
190 
191  if ((!(2&useXiLanguage)) && (!(4&useXiLanguage))) return usedV0;
192  if (useXiLanguage<6)
193  {
194  if ((2&useXiLanguage) && (v0Vertex->chiSquared()<0)) return usedV0;
195  if ((4&useXiLanguage) && (v0Vertex->chiSquared()>=0)) return usedV0;
196  }
197 
198  long myChi = -1*(long)v0Vertex->chiSquared();
199 
200  if( !useSVT && (myChi&(( long)1<<3))) return usedV0;
201  if( useSVT && !(myChi&(( long)1<<3))) return usedV0;
202 
203 
205  unsigned int k;
206  int iflag1,i,charge,tries,negKey,posKey;
207  double mlam,mala,ptV0,ptBach,radiusBach,rsq;
208  StPhysicalHelixD tmpHelix;
209  StThreeVectorF xPvx;
210  StThreeVectorD xpp,pp,impact,tmp3V,xV0,pV0,dpV0;
211  StLorentzVectorD posVec,negVec;
212 
213  //Subroutine casc_geom
214  double xc,yc,xd,yd,atmp,btmp,ctmp,dtmp,abtmp,xOut[2],yOut[2],dist;
215 
216  //Subroutine update_track_param
217  double arg,axb,ds,dz;
218  StThreeVectorF xOrig;
219 
220  //After subroutine track_mom
221  double dv0dotdb,denom,s2,valid,xAns,yAns,yy,zz,s1;
222  StThreeVectorF pBach,diffc,batv,v0atv;
223 
224  //Bloc 5bis
225  double dca,bxi;
226  double ptot_v02,bdotx,vdotx,ppar,pper;
227  StThreeVectorF pXi;
228 
229  //Function helixDCA
230  int h_tmp;
231  double pt_tmp,bcharge_tmp,curvature_tmp,dip_tmp,phase_tmp;
232  StThreeVectorD origin_tmp;
233 
234  //Rotating
235  double epsDipAngle,cstPsi;
236  StThreeVectorF epsOrigin,epsMomentum,cstOrigin;
237 
238 
239  // All the stuff without which you can't even begin to work
240  StSPtrVecXiVertex& xiVertices = *vecXi;
241  charge=0;
242  xPvx=mainv;
243  negKey=v0Vertex->daughter(negative)->key();
244  posKey=v0Vertex->daughter(positive)->key();
245  parsXi = exipar->GetTable(det_id_v0-1); //Xi cut parameters using detector id from V0
246  xV0=v0Vertex->position();
247  impact = xV0-xPvx;
248  //Cut: V0 decay length
249  if (impact.mag2() < (parsXi->rv_v0*parsXi->rv_v0)) return usedV0;
250  pV0=v0Vertex->momentum();
251  dpV0.setX(pV0.x()/abs(pV0));
252  dpV0.setY(pV0.y()/abs(pV0));
253  dpV0.setZ(pV0.z()/abs(pV0));
254  ptV0=::sqrt(pV0.x()*pV0.x()+pV0.y()*pV0.y());
255 
256  // All the stuff needed to make rotating simpler inside the loop
257  epsDipAngle=1.;
258  epsOrigin.setX(1.);
259  epsOrigin.setY(1.);
260  epsOrigin.setZ(1.);
261  epsMomentum.setX(1.);
262  epsMomentum.setY(1.);
263  epsMomentum.setZ(1.);
264  cstPsi=0.;
265  cstOrigin.setX(0.);
266  cstOrigin.setY(0.);
267  cstOrigin.setZ(0.);
268  if (1&useRotating)
269  {epsOrigin.setX(-1.);
270  epsOrigin.setY(-1.);
271  epsMomentum.setX(-1.);
272  epsMomentum.setY(-1.);
273  cstPsi=M_PI;
274  cstOrigin.setX(2*xPvx.x());
275  cstOrigin.setY(2*xPvx.y());
276  }
277  if (2&useRotating)
278  {epsDipAngle=-1.;
279  epsOrigin.setZ(-1.);
280  epsMomentum.setZ(-1.);
281  cstOrigin.setZ(2*xPvx.z());
282  }
283 
284  // Calculates Lambda invariant mass and decides if Lam or antiLam.
285  const StThreeVectorF& posVec3 = v0Vertex->momentumOfDaughter(positive);
286  const StThreeVectorF& negVec3 = v0Vertex->momentumOfDaughter(negative);
287  posVec.setVect(posVec3);
288  negVec.setVect(negVec3);
289  float pVmag2 = posVec3.mag2();
290  float nVmag2 = negVec3.mag2();
291 
292  posVec.setE(TMath::Sqrt(pVmag2+proton_mass_c2*proton_mass_c2));
293  negVec.setE(TMath::Sqrt(nVmag2+pion_minus_mass_c2*pion_minus_mass_c2));
294  mlam = (posVec+negVec).m();
295  Bool_t lamCand = (TMath::Abs(mlam-lambda_mass_c2) < parsXi->dmass);
296 
297  posVec.setE(TMath::Sqrt(pVmag2+pion_plus_mass_c2*pion_plus_mass_c2));
298  negVec.setE(TMath::Sqrt(nVmag2+proton_mass_c2*proton_mass_c2));
299  mala = (posVec+negVec).m();
300  Bool_t alaCand = (TMath::Abs(mala-lambda_mass_c2) < parsXi->dmass);
301 
302  //Cut: Lambda invariant mass
303  AssignSign:
304  if (alaCand)
305  {charge=1;
306  if (v0Vertex->dcaDaughterToPrimaryVertex(positive) < parsXi->bpn_v0)
307  {alaCand = kFALSE;
308  goto AssignSign;
309  }
310  }
311  else if (lamCand)
312  {charge=-1;
313  if (v0Vertex->dcaDaughterToPrimaryVertex(negative) < parsXi->bpn_v0)
314  {lamCand = kFALSE;
315  goto AssignSign;
316  }
317  }
318  else
319  {return usedV0;
320  }
321  charge=-(useLikesign-1)*charge;
322 
323  StHelixModel* bachGeom = new StHelixModel;
324  if (bachGeom == NULL) {gMessMgr->Info()<<"StXiFinderMaker : CAUTION : pointer bachGeom is null."<<endm; return usedV0;}
325  StHelixModel* bachGeom2 = new StHelixModel;
326  if (bachGeom2 == NULL) {gMessMgr->Info()<<"StXiFinderMaker : CAUTION : pointer bachGeom2 is null."<<endm; return usedV0;}
327 
328  // Loop over tracks (bachelors) to find Xis
329  for (k=0; k<trks; k++)
330  {bachGeom->setCharge(trk[k]->geometry()->charge());
331  bachGeom->setHelicity(trk[k]->geometry()->helicity());
332  bachGeom->setCurvature(trk[k]->geometry()->curvature());
333  bachGeom->setPsi(trk[k]->geometry()->psi()+cstPsi);
334  bachGeom->setDipAngle(epsDipAngle*trk[k]->geometry()->dipAngle());
335  bachGeom->setOrigin(cstOrigin+epsOrigin.pseudoProduct(trk[k]->geometry()->origin()));
336  bachGeom->setMomentum(epsMomentum.pseudoProduct(trk[k]->geometry()->momentum()));
337  //Check that the sign of the bachelor is the one we want.
338  if (charge*bachGeom->charge() < 0) continue;
339  //Check that ITTF and TPT tracks/V0's are not combined together.
340  if (GetTrackerUsage() == kTrackerUseBOTH)
341  {
342  //if ((v0Vertex->dcaDaughters() <= 0) && (trk[k]->fittingMethod() == TPTflag)) continue;
343  if ((v0Vertex->dcaDaughters() <= 0) && (trk[k]->fittingMethod() != ITTFflag)) continue;
344  if ((v0Vertex->dcaDaughters() >= 0) && (trk[k]->fittingMethod() == ITTFflag)) continue;
345  }
346  //Check that the bachelor is not one of the V0's daughters.
347  if (trkID[k] == negKey) continue;
348  if (trkID[k] == posKey) continue;
350 
351  //Determine detector id of pair for parsXi
352  det_id_xi=TMath::Max(det_id_v0,detId[k]);
353  //Xi cut parameters
354  parsXi=exipar->GetTable(det_id_xi-1);
355 
356  //Book a temporary clone of the bachelor helix, whose origin will be moved
357  bachGeom2->setCharge(bachGeom->charge());
358  bachGeom2->setHelicity(bachGeom->helicity());
359  bachGeom2->setCurvature(bachGeom->curvature());
360  bachGeom2->setPsi(bachGeom->psi());
361  bachGeom2->setDipAngle(bachGeom->dipAngle());
362  bachGeom2->setOrigin(bachGeom->origin());
363  bachGeom2->setMomentum(bachGeom->momentum());
364  if ((bachGeom->origin().x()==0) && (bachGeom->origin().y()==0) && (bachGeom->origin().z()==0))
365  {gMessMgr->Info()<<"StXiFinderMaker : CAUTION : bachelor candidate has all parameters = 0."<<endm;
366  continue;
367  }
368  ptBach=pt[k];
369  radiusBach=1./bachGeom->curvature();
370  rsq=radiusBach*radiusBach;
371 
372  // Calculation of the 2 intersection points between bachelor circle and V0 straight line
373  //Subroutine casc_geom
374  iflag1=2;
375  xc=bachGeom->origin().x()-bachGeom->helicity()*radiusBach*TMath::Sin(bachGeom->psi());
376  yc=bachGeom->origin().y()+bachGeom->helicity()*radiusBach*TMath::Cos(bachGeom->psi());
377  xd=xV0.x()-xc;
378  yd=xV0.y()-yc;
379  xOut[0]=0.;
380  yOut[0]=0.;
381  xOut[1]=0.;
382  yOut[1]=0.;
383  if (pV0.x() != 0)
384  {atmp=pV0.y()/pV0.x();
385  btmp=-atmp*xd+yd;
386  dtmp=atmp*atmp+1.;
387  ctmp=dtmp*rsq-btmp*btmp;
388  if (ctmp < 0)
389  {dist=fabs(xd*pV0.y()-yd*pV0.x())/ptV0;
390  if (dist >= (radiusBach+parsXi->dca_max)) iflag1=0;
391  else
392  {iflag1=1;
393  ctmp=::sqrt(rsq/(atmp*atmp+1.));
394  dtmp=-atmp*ctmp;
395  btmp=dtmp*xd+ctmp*yd;
396  if (btmp > 0.)
397  {xOut[0]=dtmp+xc;
398  yOut[0]=ctmp+yc;
399  }
400  else
401  {xOut[0]=-dtmp+xc;
402  yOut[0]=-ctmp+yc;
403  }
404  }
405  }
406  else
407  {if (ctmp == 0) iflag1=1;
408  ctmp=::sqrt(ctmp);
409  abtmp=atmp*btmp;
410  btmp=btmp+yc;
411  xOut[0]=(-abtmp+ctmp)/dtmp+xc;
412  xOut[1]=(-abtmp-ctmp)/dtmp+xc;
413  yOut[0]=atmp*(xOut[0]-xc)+btmp;
414  yOut[1]=atmp*(xOut[1]-xc)+btmp;
415  }
416  }
417  else //pV0.x()==0
418  {xOut[0]=xV0.x();
419  xOut[1]=xV0.x();
420  ctmp=rsq-xd*xd;
421  if (ctmp <= 0)
422  {dist=fabs(xd*pV0.y()-yd*pV0.x())/ptV0;
423  if (dist >= (radiusBach+parsXi->dca_max)) iflag1=0;
424  else
425  {iflag1=1;
426  yOut[0]=yc;
427  if (xV0.x() > xc) xOut[0]=xc+radiusBach;
428  else xOut[0]=xc-radiusBach;
429  }
430  }
431  else
432  {if (ctmp == 0) iflag1=1;
433  ctmp=::sqrt(ctmp);
434  yOut[0]=yc+ctmp;
435  yOut[1]=yc-ctmp;
436  }
437  }
438  //End of casc_geom
439  //Cut: drop candidate if no intersection point in 2D between V0 and bachelor trajectories
440  if (iflag1 == 0) continue; //No intersection points
441 
442  // Loop over the (1 or) 2 intersection points between bachelor circle and V0 straight line
443  for (i=0;i<iflag1;i++)
444  {//Subroutine update_track_param
445  axb=(bachGeom->origin().x()-xc)*(yOut[i]-yc)-(bachGeom->origin().y()-yc)*(xOut[i]-xc);
446  arg=axb/rsq;
447  if (arg > 1.) arg=1.;
448  if (arg < -1.) arg=-1.;
449  ds=radiusBach*TMath::ASin(arg);
450  dz=ds*TMath::Tan(bachGeom->dipAngle());
451  xOrig.setX(xOut[i]);
452  if (xOut[i] == 0.) xOrig.setX(0.01);
453  xOrig.setY(yOut[i]);
454  xOrig.setZ(bachGeom->origin().z()-(bachGeom->charge()*(Bfield/tesla)/fabs(bachGeom->charge()*(Bfield/tesla)))*dz); //Field wanted in tesla
455  bachGeom2->setOrigin(xOrig);
456  bachGeom2->setPsi(bachGeom->psi()+TMath::ASin(arg));
457  //End of update_track_param
458 
459  //Subroutine track_mom
460  xOrig=bachGeom2->momentum();
461  xOrig.setX(ptBach*TMath::Cos(bachGeom2->psi()));
462  xOrig.setY(ptBach*TMath::Sin(bachGeom2->psi()));
463  //End of track_mom
464 
465  pBach.setX(xOrig.x()/abs(xOrig));
466  pBach.setY(xOrig.y()/abs(xOrig));
467  pBach.setZ(xOrig.z()/abs(xOrig));
468  dv0dotdb=dpV0.x()*pBach.x()+dpV0.y()*pBach.y()+dpV0.z()*pBach.z();
469  diffc.setX(xV0.x()-xOut[i]);
470  diffc.setY(xV0.y()-yOut[i]);
471  diffc.setZ(xV0.z()-bachGeom2->origin().z());
472  //s1 and s2 are the distances from a point on the lines to the
473  // closest distance of approach of the lines in space
474  denom=dv0dotdb*dv0dotdb-1.;
475  s2=(dpV0.x()*dv0dotdb-pBach.x())*diffc.x() + (dpV0.y()*dv0dotdb-pBach.y())*diffc.y() + (dpV0.z()*dv0dotdb-pBach.z())*diffc.z();
476  s2=s2/denom;
477  //Check validity of linear approx. (distance moved in x and y << r1)
478  // If only mildly invalid, re-try starting with new point (up to 3 tries)
479  tries=1;
480  valid=fabs(s2*s2*(pBach.x()*pBach.x()+pBach.y()*pBach.y()));
481  //Cut: drop candidate if linear approximation too bad to determine the dca
482  while ((valid < (0.0004*rsq)) && (tries < 4) && (valid > (rsq*1.e-6)))
483  {tries++;
484  batv.setX(pBach.x()*s2+xOut[i]);
485  batv.setY(pBach.y()*s2+yOut[i]);
486  batv.setZ(pBach.z()*s2+bachGeom2->origin().z());
487 
488  //Subroutine ev0_project_track
489  dtmp=xc-batv.x();
490  atmp=yc-batv.y();
491  if (atmp == 0.)
492  {if (dtmp >= 0) xAns=xc-radiusBach;
493  else xAns=xc+radiusBach;
494  yAns=yc;
495  }
496  else
497  {ctmp=dtmp/atmp;
498  yy=radiusBach/::sqrt(ctmp*ctmp+1.);
499  zz=ctmp*yy;
500  if (atmp > 0.)
501  {xAns=-zz+xc;
502  yAns=-yy+yc;
503  }
504  else
505  {xAns=zz+xc;
506  yAns=yy+yc;
507  }
508  }
509  //End of ev0_project_track
510  xOut[i]=xAns;
511  yOut[i]=yAns;
512 
513  //Subroutine update_track_param
514  axb=(bachGeom->origin().x()-xc)*(yOut[i]-yc)-(bachGeom->origin().y()-yc)*(xOut[i]-xc);
515  arg=axb/rsq;
516  if (arg > 1.) arg=1.;
517  if (arg < -1.) arg=-1.;
518  ds=radiusBach*TMath::ASin(arg);
519  dz=ds*TMath::Tan(bachGeom->dipAngle());
520  xOrig.setX(xOut[i]);
521  if (xOut[i] == 0.) xOrig.setX(0.01);
522  xOrig.setY(yOut[i]);
523  xOrig.setZ(bachGeom->origin().z()-(bachGeom->charge()*(Bfield/tesla)/fabs(bachGeom->charge()*(Bfield/tesla)))*dz); //Field wanted in tesla
524  bachGeom2->setOrigin(xOrig);
525  bachGeom2->setPsi(bachGeom->psi()+TMath::ASin(arg));
526  //End of update_track_param
527 
528  //Subroutine track_mom
529  xOrig=bachGeom2->momentum();
530  xOrig.setX(ptBach*TMath::Cos(bachGeom2->psi()));
531  xOrig.setY(ptBach*TMath::Sin(bachGeom2->psi()));
532  //End of track_mom
533 
534  pBach.setX(xOrig.x()/abs(xOrig));
535  pBach.setY(xOrig.y()/abs(xOrig));
536  pBach.setZ(xOrig.z()/abs(xOrig));
537  dv0dotdb=dpV0.x()*pBach.x()+dpV0.y()*pBach.y()+dpV0.z()*pBach.z();
538  diffc.setX(xV0.x()-xOut[i]);
539  diffc.setY(xV0.y()-yOut[i]);
540  diffc.setZ(xV0.z()-bachGeom2->origin().z());
541  denom=dv0dotdb*dv0dotdb-1.;
542  s2=(dpV0.x()*dv0dotdb-pBach.x())*diffc.x() + (dpV0.y()*dv0dotdb-pBach.y())*diffc.y() + (dpV0.z()*dv0dotdb-pBach.z())*diffc.z();
543  s2=s2/denom;
544  valid=fabs(s2*s2*(pBach.x()*pBach.x()+pBach.y()*pBach.y()));
545  }//End of the while-loop.
546  //Cut: (same as above) drop candidate if linear approximation too bad to determine the dca
547  if ((valid < (0.0004*rsq)) && (tries < 4))
548  {batv.setX(pBach.x()*s2+xOut[i]);
549  batv.setY(pBach.y()*s2+yOut[i]);
550  batv.setZ(pBach.z()*s2+bachGeom2->origin().z());
551  //Beginning of block 5.
552  s1=(pBach.x()*dv0dotdb-dpV0.x())*diffc.x() + (pBach.y()*dv0dotdb-dpV0.y())*diffc.y() + (pBach.z()*dv0dotdb-dpV0.z())*diffc.z();
553  s1=-s1/denom;
554  v0atv.setX(dpV0.x()*s1+xV0.x());
555  v0atv.setY(dpV0.y()*s1+xV0.y());
556  v0atv.setZ(dpV0.z()*s1+xV0.z());
557  //End of block 5 and beginning of block 5bis.
558  //Cut: remove if V0 points away from Xi vertex
559  if (((xV0.x()-v0atv.x())*pV0.x() + (xV0.y()-v0atv.y())*pV0.y() + (xV0.z()-v0atv.z())*pV0.z()) <= 0.) continue;
560  dca=(v0atv.x()-batv.x())*(v0atv.x()-batv.x()) + (v0atv.y()-batv.y())*(v0atv.y()-batv.y()) + (v0atv.z()-batv.z())*(v0atv.z()-batv.z());
561  //Cut: dca Xi daughters
562  if (dca > (parsXi->dca_max*parsXi->dca_max)) continue;
563  xpp.setX((v0atv.x()+batv.x())/2.);
564  xpp.setY((v0atv.y()+batv.y())/2.);
565  xpp.setZ((v0atv.z()+batv.z())/2.);
566  //Cut: Xi decay length
567  if (((xpp.x()-xPvx.x())*(xpp.x()-xPvx.x())+(xpp.y()-xPvx.y())*(xpp.y()-xPvx.y())+(xpp.z()-xPvx.z())*(xpp.z()-xPvx.z())) <= (parsXi->rv_xi*parsXi->rv_xi)) continue;
568  //Calculate xi impact parameter
569  pXi.setX(pV0.x()+xOrig.x());
570  pXi.setY(pV0.y()+xOrig.y());
571  pXi.setZ(pV0.z()+xOrig.z());
572  //Cut: remove if Xi points away from primary vertex
573  if (((xpp.x()-xPvx.x())*pXi.x()+(xpp.y()-xPvx.y())*pXi.y()+(xpp.z()-xPvx.z())*pXi.z()) < 0.0) continue;
574  //Calculate pt-Armenteros
575  ptot_v02=ptV0*ptV0+pV0.z()*pV0.z();
576  bdotx=xOrig.x()*pXi.x()+xOrig.y()*pXi.y()+xOrig.z()*pXi.z();
577  vdotx=pV0.x()*pXi.x()+pV0.y()*pXi.y()+pV0.z()*pXi.z();
578  if (bachGeom->charge() > 0)
579  {ppar=bdotx/pXi.mag();
580  pper=(ptot[k]*ptot[k]-ppar*ppar);
581  pper=(pper>0)? ::sqrt(ptot[k]*ptot[k]-ppar*ppar):0;
582  }
583  else
584  {ppar=vdotx/pXi.mag();
585  pper=(ptot_v02-ppar*ppar);
586  pper=(pper>0)? ::sqrt(ptot_v02-ppar*ppar):0;
587  }
588  //Cut: pt-Armanteros
589  if (pper > 0.33) continue;
590  //Function helixDCA(charge,xpp,pXi,bxi);
591  //helixDCA is defined in exi_c_utils.cc (pams/global/exi/).
592  pt_tmp = ::sqrt(pXi.x()*pXi.x()+pXi.y()*pXi.y());
593  bcharge_tmp = charge*(Bfield/kilogauss); //Field wanted in kGauss
594  curvature_tmp = TMath::Abs(bcharge_tmp)*C_D_CURVATURE/pt_tmp;
595  dip_tmp = atan(pXi.z()/pt_tmp);
596  h_tmp = ((bcharge_tmp > 0) ? -1 : 1);
597  phase_tmp = atan2(pXi.y(),pXi.x())-(h_tmp*M_PI_2);
598  origin_tmp.setX(xpp.x());
599  origin_tmp.setY(xpp.y());
600  origin_tmp.setZ(xpp.z());
601  StHelixD *globHelix = new StHelixD(curvature_tmp, dip_tmp, phase_tmp, origin_tmp, h_tmp);
602  bxi=globHelix->distance(xPvx);
603  //End of helixDCA
604  //Cut: dca Xi to primary vertex
605  if (bxi > parsXi->bxi_max)
606  {delete globHelix;
607  globHelix=0;
608  continue;
609  }
610  StThreeVectorD p1 = globHelix->at(globHelix->pathLength(xPvx));
611  StThreeVectorD p2(p1.x()-globHelix->xcenter(),p1.y()-globHelix->ycenter(),0);
612  StThreeVectorD p3(xPvx.x()-globHelix->xcenter(),xPvx.y()-globHelix->ycenter(),0);
613  if (p3.mag2() > p2.mag2()) bxi=-bxi;
614  delete globHelix;
615  globHelix=0;
616  xiVertex = new StXiVertex();
617  xiVertex->setPosition(xpp);
618  xiVertex->addDaughter(trk[k]);
619  xiVertex->setDcaBachelorToPrimaryVertex(trk[k]->impactParameter());
620  xiVertex->setMomentumOfBachelor(xOrig);
621  xiVertex->setDcaDaughters(::sqrt(dca));
622  xiVertex->setDcaParentToPrimaryVertex(bxi);
623  xiVertex->setV0Vertex(v0Vertex);
624 
626  //Set chi2 to trace SVT usage
627  long int v0ChiSq = -1*(long int)v0Vertex->chiSquared();
628  if(detId[k]==2 || detId[k]==3)
629  {//if an SVT track was used on bachelor set the last bit to 1
630  v0ChiSq |=((long int)1 << 0);
631  }
632  v0ChiSq *=-1;
633  xiVertex->setChiSquared((float)v0ChiSq);
635 
636  xiVertices.push_back(xiVertex);
637  usedV0 = kTRUE;
638  } //End if (valid && tries<4)
639  } //End loop over 2 intersection points
640  } // k-Loop
641  delete bachGeom;
642  bachGeom=0;
643  delete bachGeom2;
644  bachGeom2=0;
645 
646  if (alaCand)
647  {alaCand = kFALSE;
648  goto AssignSign;
649  }
650 
651  return usedV0;
652 }
653 //_____________________________________________________________________________
654 // $Id: StXiFinderMaker.cxx,v 1.25 2017/01/06 21:01:50 smirnovd Exp $
655 // $Log: StXiFinderMaker.cxx,v $
656 // Revision 1.25 2017/01/06 21:01:50 smirnovd
657 // Use pi constant from standard library, s/C_PI/M_PI/
658 //
659 // Revision 1.24 2016/12/12 17:18:04 smirnovd
660 // Removed outdated ClassImp ROOT macro
661 //
662 // Revision 1.23 2008/04/03 19:58:36 fisyak
663 // move parameters initialization from Init into InitRun
664 //
665 // Revision 1.22 2005/07/19 22:10:14 perev
666 // STARFPE
667 //
668 // Revision 1.21 2004/04/15 20:15:53 jeromel
669 // Forgot one of them
670 //
671 // Revision 1.20 2004/04/02 08:57:34 faivre
672 // Use actual TPT flag rather than "not ITTF" for TPT tracks.
673 //
674 // Revision 1.19 2004/02/05 16:52:28 faivre
675 // Add cut pt-Armanteros > 0.33 ; small cleanup.
676 //
677 // Revision 1.18 2004/02/04 14:53:36 faivre
678 // Remove bug introduced in previous version (confusion dca/bxi when copy-pasting from helixDCA).
679 //
680 // Revision 1.17 2004/02/04 14:24:58 faivre
681 // Add sign of dcaXiDaughters, slightly move cuts, small cleanup.
682 //
683 // Revision 1.16 2004/02/03 14:29:36 faivre
684 // Spring-cleaning 2 months early ;-) Algo strictly equivalent to previous one although huge reshaping.
685 //
686 // Revision 1.15 2004/02/02 12:10:54 faivre
687 // XiFinder now able to run on muDsts as well :-) + update user-friendliness.
688 //
689 // Revision 1.14 2003/11/08 18:26:00 faivre
690 // Bfield + consistency int/short
691 //
692 // Revision 1.13 2003/09/17 12:00:23 faivre
693 // RH8 : initialize everything in constructor.
694 //
695 // Revision 1.12 2003/09/02 17:58:59 perev
696 // gcc 3.2 updates + WarnOff
697 //
698 // Revision 1.11 2003/08/22 17:47:14 caines
699 // Get sign AND magnitude of mag field correctly for Xi and V0 finder
700 //
701 // Revision 1.10 2003/07/15 17:39:16 faivre
702 // Doesn't take Bfield from gufld anymore (takes Bfield calculated in V0Finder).
703 //
704 // Revision 1.9 2003/07/04 17:52:46 faivre
705 // Use SVT cuts if any dg has a SVT hit.
706 //
707 // Revision 1.8 2003/06/24 16:20:11 faivre
708 // Uses SVT tracks. Fixed bool calculations. Exits when bad param. Reshaping.
709 //
710 // Revision 1.7 2003/05/14 19:15:36 faivre
711 // Fancy choices Fortran/C++ V0's and Xi's. Xi rotating and like-sign.
712 //
713 // Revision 1.6 2003/05/07 10:51:42 faivre
714 // Use brand new StHelixModel::setMomentum to solve memory leaks.
715 //
716 // Revision 1.4 2003/05/02 21:21:08 lbarnby
717 // Now identify ITTF tracks by fittingMethod() equal to kITKalmanFitId
718 //
719 // Revision 1.3 2003/04/30 20:38:22 perev
720 // Warnings cleanup. Modified lines marked VP
721 //
722 // Revision 1.2 2003/04/30 19:16:04 faivre
723 // Fix storage part. ITTF vs TPT Xis.
724 //
725 //
StV0FinderMaker::Prepare
virtual Int_t Prepare()
Definition: StV0FinderMaker.cxx:203
StXiFinderMaker::UseV0
virtual Bool_t UseV0()
Definition: StXiFinderMaker.cxx:187
StXiFinderMaker::Make
virtual Int_t Make()
Definition: StXiFinderMaker.cxx:135
StHelix::distance
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
StV0FinderMaker::Make
virtual Int_t Make()
Definition: StV0FinderMaker.cxx:392
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StHelix::ycenter
double ycenter() const
x-center of circle in xy-plane
Definition: StHelix.cc:215
kStErr
Definition: Stypes.h:44
StHelix::xcenter
double xcenter() const
aziumth in xy-plane measured from ring center
Definition: StHelix.cc:207
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5