StRoot: StSecondaryVertexMaker/StXiFinderMaker.cxx Source File

00001 
00002 
00003 
00004 
00005 
00006 //
00007 //  Cuts can be found in the code by comments beginning with "Cut:"
00008 //
00009 //
00010 
00011 
00012 #include "StXiFinderMaker.h"
00013 #include "StMessMgr.h"
00014 #include "StEvent/StEventTypes.h"
00015 #include "TMath.h"
00016 #include "TVector2.h"
00017 #include "tables/St_exi_exipar_Table.h"
00018 #include "PhysicalConstants.h"
00019 #include "math_constants.h"
00020 #include "phys_constants.h"
00021 
00022 
00023 StSPtrVecXiVertex* vecXi=0;
00024 
00025 ClassImp(StXiFinderMaker)
00026 
00027 
00028 //_____________________________________________________________________________
00029 StXiFinderMaker::StXiFinderMaker(const char *name):StV0FinderMaker(name),
00030 exipar(0),parsXi(0),xiVertex(0),det_id_xi(0){}
00031 
00032 //_____________________________________________________________________________
00033 StXiFinderMaker::~StXiFinderMaker() {}
00034 //_____________________________________________________________________________
00035 Int_t StXiFinderMaker::Init()
00036 {bool a,b,c;
00037  if ((useTracker!=kTrackerUseTPT) && (useTracker!=kTrackerUseITTF) && (useTracker!=kTrackerUseBOTH))
00038     {gMessMgr->Error("StXiFinderMaker::Init() : wrong TrackerUsage parameter set.");
00039      return kStErr;
00040      }
00041  if ((useSVT!=kNoSVT) && (useSVT!=kUseSVT))
00042     {gMessMgr->Error("StXiFinderMaker::Init() : wrong SVTUsage parameter set.");
00043      return kStErr;
00044      }
00045  if ((useEventModel!=kUseStEvent) && (useEventModel!=kUseMuDst))
00046     {gMessMgr->Error("StXiFinderMaker::Init() : wrong EventModelUsage parameter set.");
00047      return kStErr;
00048      }
00049  if ((useLikesign!=kLikesignUseStandard) && (useLikesign!=kLikesignUseLikesign))
00050     {gMessMgr->Error("StXiFinderMaker::Init() : wrong LikesignUsage parameter set.");
00051      return kStErr;
00052      }
00053  if ((useRotating!=kRotatingUseStandard) && (useRotating!=kRotatingUseRotating) && (useRotating!=kRotatingUseSymmetry) && (useRotating!=kRotatingUseRotatingAndSymmetry))
00054     {gMessMgr->Error("StXiFinderMaker::Init() : wrong RotatingUsage parameter set.");
00055      return kStErr;
00056      }
00057  
00058  if (useTracker == kTrackerUseTPT) gMessMgr->Info()<<"StXiFinderMaker : use TPT tracks."<<endm;
00059  if (useTracker == kTrackerUseITTF) gMessMgr->Info()<<"StXiFinderMaker : use ITTF tracks."<<endm;
00060  if (useTracker == kTrackerUseBOTH) gMessMgr->Info()<<"StXiFinderMaker : use TPT *and* ITTF tracks."<<endm;
00061  if (useSVT == kUseSVT) gMessMgr->Info()<<"StXiFinderMaker : use SVT points if possible."<<endm;
00062  if (useSVT == kNoSVT) gMessMgr->Info()<<"StXiFinderMaker : do not use SVT points."<<endm;
00063  if (useEventModel == kUseStEvent) gMessMgr->Info()<<"StXiFinderMaker : expect StEvent files in input."<<endm;
00064  if (useEventModel == kUseMuDst)  gMessMgr->Info()<<"StXiFinderMaker : expect MuDst files in input."<<endm;
00065  if (useLikesign == kLikesignUseLikesign) gMessMgr->Info()<<"StXiFinderMaker : does like-sign finding."<<endm;
00066  if (useRotating == kRotatingUseRotating) gMessMgr->Info()<<"StXiFinderMaker : does rotating finding."<<endm;
00067  if (useRotating == kRotatingUseSymmetry) gMessMgr->Info()<<"StXiFinderMaker : does symmetry finding."<<endm;
00068  if (useRotating == kRotatingUseRotatingAndSymmetry) gMessMgr->Info()<<"StXiFinderMaker : does rotating + symmetry finding."<<endm;
00069  
00070  if (useLanguage != kLanguageUseSpecial)
00071     {a=(bool)(1&(useLanguage>>2));
00072      b=(bool)(1&(useLanguage>>1));
00073      c=(bool)(1&useLanguage);
00074      useV0Language=2*(!(a^c))+(a|c);
00075      useXiLanguage=4*(b&(!(a^c)))+2*(a&b&(!c))+(a|c);
00076      }
00077  switch (useLanguage)
00078     {case kLanguageUseOldRun : gMessMgr->Info()<<"StXiFinderMaker : Fortran run."<<endm;
00079                                break;
00080      case kLanguageUseRun : gMessMgr->Info()<<"StXiFinderMaker : C++ run."<<endm;
00081                             break;
00082      case kLanguageUseTestV0Finder : gMessMgr->Info()<<"StXiFinderMaker : Test V0Finder."<<endm;
00083                                      break;
00084      case kLanguageUseTestXiFinder : gMessMgr->Info()<<"StXiFinderMaker : Test XiFinder."<<endm;
00085                                      break;
00086      case kLanguageUseTestBothFinders : gMessMgr->Info()<<"StXiFinderMaker : Test V0Finder and XiFinder."<<endm;
00087                                         break;
00088      case kLanguageUseSpecial : break;
00089      default : gMessMgr->Error("StXiFinderMaker::Init() : wrong LanguageUsage parameter set.");
00090                return kStErr;
00091      }
00092  if ((useXiLanguage!=kXiLanguageUseFortran) && (useXiLanguage!=kXiLanguageUseCppOnFortranV0) && (useXiLanguage!=kXiLanguageUseCppOnCppV0) && (useXiLanguage!=kXiLanguageUseFortranAndCppOnFortranV0) && (useXiLanguage!=kXiLanguageUseFortranAndCppOnCppV0) && (useXiLanguage!=kXiLanguageUseBothCpp) && (useXiLanguage!=kXiLanguageUseAll))
00093     {gMessMgr->Error("StXiFinderMaker::Init() : wrong XiLanguageUsage parameter set.");
00094      return kStErr;
00095      }
00096  switch (useV0Language)
00097     {case kV0LanguageUseFortran : if ((useXiLanguage!=kXiLanguageUseFortran) && (useXiLanguage!=kXiLanguageUseCppOnFortranV0) && (useXiLanguage!=kXiLanguageUseFortranAndCppOnFortranV0))
00098                                      {gMessMgr->Info()<<"StXiFinderMaker : BE CAREFUL : impossible combination asked."<<endm;
00099                                       gMessMgr->Info()<<"StXiFinderMaker :    Set it to testXiFinder."<<endm;
00100                                       useXiLanguage=kXiLanguageUseFortranAndCppOnFortranV0;
00101                                       }
00102                                   break;
00103      case kV0LanguageUseCpp : if (useXiLanguage!=kXiLanguageUseCppOnCppV0)
00104                                  {gMessMgr->Info()<<"StXiFinderMaker : BE CAREFUL : impossible combination asked."<<endm;
00105                                   gMessMgr->Info()<<"StXiFinderMaker :    Set it to normalRun."<<endm;
00106                                   useXiLanguage=kXiLanguageUseCppOnCppV0;
00107                                   }
00108                               break;
00109      case kV0LanguageUseBoth : break;
00110      default : gMessMgr->Error("StXiFinderMaker::Init() : wrong V0LanguageUsage parameter set.");
00111                return kStErr;
00112      }
00113  if (1&useV0Language) gMessMgr->Info()<<"StXiFinderMaker :    Will store Fortran V0s."<<endm;
00114  if (2&useV0Language) gMessMgr->Info()<<"StXiFinderMaker :    Will store C++ V0s."<<endm;
00115  if (1&useXiLanguage) gMessMgr->Info()<<"StXiFinderMaker :    Will store Fortran Xis."<<endm;
00116  if (2&useXiLanguage) gMessMgr->Info()<<"StXiFinderMaker :    Will store C++ Xis made with Fortran V0s."<<endm;
00117  if (4&useXiLanguage) gMessMgr->Info()<<"StXiFinderMaker :    Will store C++ Xis made with C++ V0s."<<endm;
00118 
00119  if (useEventModel) //initialize mMuDstMaker
00120     {mMuDstMaker = (StMuDstMaker*)GetMaker("myMuDstMaker");
00121      if(!mMuDstMaker) gMessMgr->Warning("StXiFinderMaker::Init can't find a valid MuDst.");
00122      }
00123  
00124  return StMaker::Init();
00125  }
00126 //=============================================================================
00127 Int_t StXiFinderMaker::InitRun(int runumber) {
00128  exipar = (St_exi_exipar*) GetDataBase("global/vertices/exipar");
00129  if (!exipar)
00130     {gMessMgr->Error("StXiFinderMaker::Init() : could not find exipar in database.");
00131      return kStErr;
00132      }
00133   return StMaker::InitRun(runumber);
00134  }
00135 //=============================================================================
00136 //_____________________________________________________________________________
00137 Int_t StXiFinderMaker::Make() {
00138 
00139   Int_t iRes;
00140 
00141   // Prepare event and track variables
00142   iRes = Prepare();
00143   if (iRes != kStOk) return iRes;
00144 
00145   StSPtrVecXiVertex& xiVertices = event->xiVertices();
00146   gMessMgr->Info()<<"StXiFinderMaker : coming in I have "<<xiVertices.size()<<" Xis."<<endm;
00159   if (!(1&useXiLanguage) && !((4&useXiLanguage) && (!(1&useV0Language))))
00160      {// Erase existing Xis
00161       gMessMgr->Info() << "StXiFinderMaker : pre-existing Xis deleted." << endm;
00162       StSPtrVecXiVertex xiVertices2;
00163       xiVertices = xiVertices2;
00164       }
00165   vecXi = &xiVertices;
00166 
00167   // Call the V0-finding, which will in turn call
00168   // the UseV0() member function for each V0 found
00169   
00170   if (2&useXiLanguage)
00171      {StSPtrVecV0Vertex& v0Vertices = event->v0Vertices();
00172      unsigned int nV0s = v0Vertices.size();
00173      det_id_v0 = 1; // for lack of any further information
00174      for (unsigned int i=0; i<nV0s; i++)
00175         {v0Vertex = v0Vertices[i];
00176          if (v0Vertex) UseV0();
00177          }
00178      }
00179   if ((4&useXiLanguage) || (2&useV0Language))
00180     {iRes = StV0FinderMaker::Make();
00181      if (iRes != kStOk) return iRes;
00182      }
00183 
00184   gMessMgr->Info()<<"StXiFinderMaker : now I have "<<xiVertices.size()<<" Xis."<<endm;
00185 
00186   return kStOk;
00187 }
00188 //_____________________________________________________________________________
00189 Bool_t StXiFinderMaker::UseV0() {
00190   
00191   Bool_t usedV0 = kFALSE;
00192   
00193   if ((!(2&useXiLanguage)) && (!(4&useXiLanguage))) return usedV0;
00194   if (useXiLanguage<6)
00195     {
00196       if ((2&useXiLanguage) && (v0Vertex->chiSquared()<0)) return usedV0;
00197       if ((4&useXiLanguage) && (v0Vertex->chiSquared()>=0)) return usedV0;
00198     }
00199   
00200   long myChi = -1*(long)v0Vertex->chiSquared();
00201   
00202   if( !useSVT && (myChi&(( long)1<<3))) return usedV0;
00203   if( useSVT && !(myChi&(( long)1<<3))) return usedV0;
00204 
00205 
00207   unsigned int k;
00208   int iflag1,i,charge,tries,negKey,posKey;
00209   double mlam,mala,ptV0,ptBach,radiusBach,rsq;
00210   StPhysicalHelixD tmpHelix;
00211   StThreeVectorF xPvx;
00212   StThreeVectorD xpp,pp,impact,tmp3V,xV0,pV0,dpV0;
00213   StLorentzVectorD posVec,negVec;
00214   
00215   //Subroutine casc_geom
00216   double xc,yc,xd,yd,atmp,btmp,ctmp,dtmp,abtmp,xOut[2],yOut[2],dist;
00217   
00218   //Subroutine update_track_param
00219   double arg,axb,ds,dz;
00220   StThreeVectorF xOrig;
00221   
00222   //After subroutine track_mom
00223   double dv0dotdb,denom,s2,valid,xAns,yAns,yy,zz,s1;
00224   StThreeVectorF pBach,diffc,batv,v0atv;
00225   
00226   //Bloc 5bis
00227   double dca,bxi;
00228   double ptot_v02,bdotx,vdotx,ppar,pper;
00229   StThreeVectorF pXi;
00230   
00231   //Function helixDCA
00232   int h_tmp;
00233   double pt_tmp,bcharge_tmp,curvature_tmp,dip_tmp,phase_tmp;
00234   StThreeVectorD origin_tmp;
00235   
00236   //Rotating
00237   double epsDipAngle,cstPsi;
00238   StThreeVectorF epsOrigin,epsMomentum,cstOrigin;
00239   
00240 
00241   // All the stuff without which you can't even begin to work
00242   StSPtrVecXiVertex& xiVertices = *vecXi;
00243   charge=0;
00244   xPvx=mainv;
00245   negKey=v0Vertex->daughter(negative)->key();
00246   posKey=v0Vertex->daughter(positive)->key();
00247   parsXi = exipar->GetTable(det_id_v0-1); //Xi cut parameters using detector id from V0
00248   xV0=v0Vertex->position();
00249   impact = xV0-xPvx;
00250   //Cut: V0 decay length
00251   if (impact.mag2() < (parsXi->rv_v0*parsXi->rv_v0)) return usedV0;
00252   pV0=v0Vertex->momentum();
00253   dpV0.setX(pV0.x()/abs(pV0));
00254   dpV0.setY(pV0.y()/abs(pV0));
00255   dpV0.setZ(pV0.z()/abs(pV0));
00256   ptV0=::sqrt(pV0.x()*pV0.x()+pV0.y()*pV0.y());
00257   
00258   // All the stuff needed to make rotating simpler inside the loop
00259   epsDipAngle=1.;
00260   epsOrigin.setX(1.);
00261   epsOrigin.setY(1.);
00262   epsOrigin.setZ(1.);
00263   epsMomentum.setX(1.);
00264   epsMomentum.setY(1.);
00265   epsMomentum.setZ(1.);
00266   cstPsi=0.;
00267   cstOrigin.setX(0.);
00268   cstOrigin.setY(0.);
00269   cstOrigin.setZ(0.);
00270   if (1&useRotating)
00271      {epsOrigin.setX(-1.);
00272       epsOrigin.setY(-1.);
00273       epsMomentum.setX(-1.);
00274       epsMomentum.setY(-1.);
00275       cstPsi=C_PI;
00276       cstOrigin.setX(2*xPvx.x());
00277       cstOrigin.setY(2*xPvx.y());
00278       }
00279   if (2&useRotating)
00280      {epsDipAngle=-1.;
00281       epsOrigin.setZ(-1.);
00282       epsMomentum.setZ(-1.);
00283       cstOrigin.setZ(2*xPvx.z());
00284       }
00285 
00286   // Calculates Lambda invariant mass and decides if Lam or antiLam.
00287   const StThreeVectorF& posVec3 = v0Vertex->momentumOfDaughter(positive);
00288   const StThreeVectorF& negVec3 = v0Vertex->momentumOfDaughter(negative);
00289   posVec.setVect(posVec3);
00290   negVec.setVect(negVec3);
00291   float pVmag2 = posVec3.mag2();
00292   float nVmag2 = negVec3.mag2();
00293 
00294   posVec.setE(TMath::Sqrt(pVmag2+proton_mass_c2*proton_mass_c2));
00295   negVec.setE(TMath::Sqrt(nVmag2+pion_minus_mass_c2*pion_minus_mass_c2));
00296   mlam = (posVec+negVec).m();
00297   Bool_t lamCand = (TMath::Abs(mlam-lambda_mass_c2) < parsXi->dmass);
00298 
00299   posVec.setE(TMath::Sqrt(pVmag2+pion_plus_mass_c2*pion_plus_mass_c2));
00300   negVec.setE(TMath::Sqrt(nVmag2+proton_mass_c2*proton_mass_c2));
00301   mala = (posVec+negVec).m();
00302   Bool_t alaCand = (TMath::Abs(mala-lambda_mass_c2) < parsXi->dmass);
00303   
00304   //Cut: Lambda invariant mass
00305   AssignSign:
00306   if (alaCand)
00307      {charge=1;
00308       if (v0Vertex->dcaDaughterToPrimaryVertex(positive) < parsXi->bpn_v0) 
00309          {alaCand = kFALSE;
00310           goto AssignSign;
00311           }
00312       }
00313       else if (lamCand)
00314      {charge=-1;
00315       if (v0Vertex->dcaDaughterToPrimaryVertex(negative) < parsXi->bpn_v0)
00316          {lamCand = kFALSE;
00317           goto AssignSign;
00318           }
00319       }
00320       else
00321      {return usedV0;
00322       }
00323   charge=-(useLikesign-1)*charge;
00324 
00325   StHelixModel* bachGeom = new StHelixModel;
00326   if (bachGeom == NULL) {gMessMgr->Info()<<"StXiFinderMaker : CAUTION : pointer bachGeom is null."<<endm; return usedV0;}
00327   StHelixModel* bachGeom2 = new StHelixModel;
00328   if (bachGeom2 == NULL) {gMessMgr->Info()<<"StXiFinderMaker : CAUTION : pointer bachGeom2 is null."<<endm; return usedV0;}
00329   
00330   // Loop over tracks (bachelors) to find Xis
00331   for (k=0; k<trks; k++)
00332      {bachGeom->setCharge(trk[k]->geometry()->charge());
00333       bachGeom->setHelicity(trk[k]->geometry()->helicity());
00334       bachGeom->setCurvature(trk[k]->geometry()->curvature());
00335       bachGeom->setPsi(trk[k]->geometry()->psi()+cstPsi);
00336       bachGeom->setDipAngle(epsDipAngle*trk[k]->geometry()->dipAngle());
00337       bachGeom->setOrigin(cstOrigin+epsOrigin.pseudoProduct(trk[k]->geometry()->origin()));
00338       bachGeom->setMomentum(epsMomentum.pseudoProduct(trk[k]->geometry()->momentum()));
00339       //Check that the sign of the bachelor is the one we want.
00340       if (charge*bachGeom->charge() < 0) continue;
00341       //Check that ITTF and TPT tracks/V0's are not combined together.
00342       if (GetTrackerUsage() == kTrackerUseBOTH)
00343          {
00344            //if ((v0Vertex->dcaDaughters() <= 0) && (trk[k]->fittingMethod() == TPTflag)) continue;
00345            if ((v0Vertex->dcaDaughters() <= 0) && (trk[k]->fittingMethod() != ITTFflag)) continue;
00346            if ((v0Vertex->dcaDaughters() >= 0) && (trk[k]->fittingMethod() == ITTFflag)) continue;
00347          }
00348       //Check that the bachelor is not one of the V0's daughters.
00349       if (trkID[k] == negKey) continue;
00350       if (trkID[k] == posKey) continue;
00352       
00353       //Determine detector id of pair for parsXi
00354       det_id_xi=TMath::Max(det_id_v0,detId[k]);
00355       //Xi cut parameters
00356       parsXi=exipar->GetTable(det_id_xi-1);
00357       
00358       //Book a temporary clone of the bachelor helix, whose origin will be moved
00359       bachGeom2->setCharge(bachGeom->charge());
00360       bachGeom2->setHelicity(bachGeom->helicity());
00361       bachGeom2->setCurvature(bachGeom->curvature());
00362       bachGeom2->setPsi(bachGeom->psi());
00363       bachGeom2->setDipAngle(bachGeom->dipAngle());
00364       bachGeom2->setOrigin(bachGeom->origin());
00365       bachGeom2->setMomentum(bachGeom->momentum());
00366       if ((bachGeom->origin().x()==0) && (bachGeom->origin().y()==0) && (bachGeom->origin().z()==0))
00367          {gMessMgr->Info()<<"StXiFinderMaker : CAUTION : bachelor candidate has all parameters = 0."<<endm;
00368           continue;
00369           }
00370       ptBach=pt[k];
00371       radiusBach=1./bachGeom->curvature();
00372       rsq=radiusBach*radiusBach;
00373       
00374       // Calculation of the 2 intersection points between bachelor circle and V0 straight line
00375       //Subroutine casc_geom
00376       iflag1=2;
00377       xc=bachGeom->origin().x()-bachGeom->helicity()*radiusBach*TMath::Sin(bachGeom->psi());
00378       yc=bachGeom->origin().y()+bachGeom->helicity()*radiusBach*TMath::Cos(bachGeom->psi());
00379       xd=xV0.x()-xc;
00380       yd=xV0.y()-yc;
00381       xOut[0]=0.;
00382       yOut[0]=0.;
00383       xOut[1]=0.;
00384       yOut[1]=0.;
00385       if (pV0.x() != 0)
00386          {atmp=pV0.y()/pV0.x();
00387           btmp=-atmp*xd+yd;
00388           dtmp=atmp*atmp+1.;
00389           ctmp=dtmp*rsq-btmp*btmp;
00390           if (ctmp < 0)
00391              {dist=fabs(xd*pV0.y()-yd*pV0.x())/ptV0;
00392               if (dist >= (radiusBach+parsXi->dca_max)) iflag1=0;
00393                   else
00394                  {iflag1=1;
00395                   ctmp=::sqrt(rsq/(atmp*atmp+1.));
00396                   dtmp=-atmp*ctmp;
00397                   btmp=dtmp*xd+ctmp*yd;
00398                   if (btmp > 0.)
00399                      {xOut[0]=dtmp+xc;
00400                       yOut[0]=ctmp+yc;
00401                       }
00402                       else
00403                      {xOut[0]=-dtmp+xc;
00404                       yOut[0]=-ctmp+yc;
00405                       }
00406                   }
00407               }
00408               else
00409              {if (ctmp == 0) iflag1=1;
00410               ctmp=::sqrt(ctmp);
00411               abtmp=atmp*btmp;
00412               btmp=btmp+yc;
00413               xOut[0]=(-abtmp+ctmp)/dtmp+xc;
00414               xOut[1]=(-abtmp-ctmp)/dtmp+xc;
00415               yOut[0]=atmp*(xOut[0]-xc)+btmp;
00416               yOut[1]=atmp*(xOut[1]-xc)+btmp;
00417               }
00418           }
00419           else //pV0.x()==0
00420          {xOut[0]=xV0.x();
00421           xOut[1]=xV0.x();
00422           ctmp=rsq-xd*xd;
00423           if (ctmp <= 0)
00424              {dist=fabs(xd*pV0.y()-yd*pV0.x())/ptV0;
00425               if (dist >= (radiusBach+parsXi->dca_max)) iflag1=0;
00426                   else
00427                  {iflag1=1;
00428                   yOut[0]=yc;
00429                   if (xV0.x() > xc) xOut[0]=xc+radiusBach;
00430                      else xOut[0]=xc-radiusBach;
00431                   }
00432               }
00433               else
00434              {if (ctmp == 0) iflag1=1;
00435               ctmp=::sqrt(ctmp);
00436               yOut[0]=yc+ctmp;
00437               yOut[1]=yc-ctmp;
00438               }
00439           }
00440       //End of casc_geom
00441       //Cut: drop candidate if no intersection point in 2D between V0 and bachelor trajectories
00442       if (iflag1 == 0) continue; //No intersection points
00443       
00444       // Loop over the (1 or) 2 intersection points between bachelor circle and V0 straight line
00445       for (i=0;i<iflag1;i++)
00446          {//Subroutine update_track_param
00447           axb=(bachGeom->origin().x()-xc)*(yOut[i]-yc)-(bachGeom->origin().y()-yc)*(xOut[i]-xc);
00448           arg=axb/rsq;
00449           if (arg > 1.) arg=1.;
00450           if (arg < -1.) arg=-1.;
00451           ds=radiusBach*TMath::ASin(arg);
00452           dz=ds*TMath::Tan(bachGeom->dipAngle());
00453           xOrig.setX(xOut[i]);
00454           if (xOut[i] == 0.) xOrig.setX(0.01);
00455           xOrig.setY(yOut[i]);
00456           xOrig.setZ(bachGeom->origin().z()-(bachGeom->charge()*(Bfield/tesla)/fabs(bachGeom->charge()*(Bfield/tesla)))*dz); //Field wanted in tesla
00457           bachGeom2->setOrigin(xOrig);
00458           bachGeom2->setPsi(bachGeom->psi()+TMath::ASin(arg));
00459           //End of update_track_param
00460           
00461           //Subroutine track_mom
00462           xOrig=bachGeom2->momentum();
00463           xOrig.setX(ptBach*TMath::Cos(bachGeom2->psi()));
00464           xOrig.setY(ptBach*TMath::Sin(bachGeom2->psi()));
00465           //End of track_mom
00466           
00467           pBach.setX(xOrig.x()/abs(xOrig));
00468           pBach.setY(xOrig.y()/abs(xOrig));
00469           pBach.setZ(xOrig.z()/abs(xOrig));
00470           dv0dotdb=dpV0.x()*pBach.x()+dpV0.y()*pBach.y()+dpV0.z()*pBach.z();
00471           diffc.setX(xV0.x()-xOut[i]);
00472           diffc.setY(xV0.y()-yOut[i]);
00473           diffc.setZ(xV0.z()-bachGeom2->origin().z());
00474           //s1 and s2 are the distances from a point on the lines to the
00475           // closest distance of approach of the lines in space
00476           denom=dv0dotdb*dv0dotdb-1.;
00477           s2=(dpV0.x()*dv0dotdb-pBach.x())*diffc.x() + (dpV0.y()*dv0dotdb-pBach.y())*diffc.y() + (dpV0.z()*dv0dotdb-pBach.z())*diffc.z();
00478           s2=s2/denom;
00479           //Check validity of linear approx. (distance moved in x and y << r1)
00480           // If only mildly invalid, re-try starting with new point (up to 3 tries)
00481           tries=1;
00482           valid=fabs(s2*s2*(pBach.x()*pBach.x()+pBach.y()*pBach.y()));
00483           //Cut: drop candidate if linear approximation too bad to determine the dca
00484           while ((valid < (0.0004*rsq)) && (tries < 4) && (valid > (rsq*1.e-6)))
00485              {tries++;
00486               batv.setX(pBach.x()*s2+xOut[i]);
00487               batv.setY(pBach.y()*s2+yOut[i]);
00488               batv.setZ(pBach.z()*s2+bachGeom2->origin().z());
00489               
00490               //Subroutine ev0_project_track
00491               dtmp=xc-batv.x();
00492               atmp=yc-batv.y();
00493               if (atmp == 0.)
00494                  {if (dtmp >= 0) xAns=xc-radiusBach;
00495                      else xAns=xc+radiusBach;
00496                   yAns=yc;
00497                   }
00498                   else
00499                  {ctmp=dtmp/atmp;
00500                   yy=radiusBach/::sqrt(ctmp*ctmp+1.);
00501                   zz=ctmp*yy;
00502                   if (atmp > 0.)
00503                      {xAns=-zz+xc;
00504                       yAns=-yy+yc;
00505                       }
00506                       else
00507                      {xAns=zz+xc;
00508                       yAns=yy+yc;
00509                       }
00510                   }
00511               //End of ev0_project_track
00512               xOut[i]=xAns;
00513               yOut[i]=yAns;
00514               
00515               //Subroutine update_track_param
00516               axb=(bachGeom->origin().x()-xc)*(yOut[i]-yc)-(bachGeom->origin().y()-yc)*(xOut[i]-xc);
00517               arg=axb/rsq;
00518               if (arg > 1.) arg=1.;
00519               if (arg < -1.) arg=-1.;
00520               ds=radiusBach*TMath::ASin(arg);
00521               dz=ds*TMath::Tan(bachGeom->dipAngle());
00522               xOrig.setX(xOut[i]);
00523               if (xOut[i] == 0.) xOrig.setX(0.01);
00524               xOrig.setY(yOut[i]);
00525               xOrig.setZ(bachGeom->origin().z()-(bachGeom->charge()*(Bfield/tesla)/fabs(bachGeom->charge()*(Bfield/tesla)))*dz); //Field wanted in tesla
00526               bachGeom2->setOrigin(xOrig);
00527               bachGeom2->setPsi(bachGeom->psi()+TMath::ASin(arg));
00528               //End of update_track_param                  
00529 
00530               //Subroutine track_mom
00531               xOrig=bachGeom2->momentum();
00532               xOrig.setX(ptBach*TMath::Cos(bachGeom2->psi()));
00533               xOrig.setY(ptBach*TMath::Sin(bachGeom2->psi()));
00534               //End of track_mom
00535               
00536               pBach.setX(xOrig.x()/abs(xOrig));
00537               pBach.setY(xOrig.y()/abs(xOrig));
00538               pBach.setZ(xOrig.z()/abs(xOrig));
00539               dv0dotdb=dpV0.x()*pBach.x()+dpV0.y()*pBach.y()+dpV0.z()*pBach.z();
00540               diffc.setX(xV0.x()-xOut[i]);
00541               diffc.setY(xV0.y()-yOut[i]);
00542               diffc.setZ(xV0.z()-bachGeom2->origin().z());
00543               denom=dv0dotdb*dv0dotdb-1.;
00544               s2=(dpV0.x()*dv0dotdb-pBach.x())*diffc.x() + (dpV0.y()*dv0dotdb-pBach.y())*diffc.y() + (dpV0.z()*dv0dotdb-pBach.z())*diffc.z();
00545               s2=s2/denom;
00546               valid=fabs(s2*s2*(pBach.x()*pBach.x()+pBach.y()*pBach.y()));
00547               }//End of the while-loop.
00548           //Cut: (same as above) drop candidate if linear approximation too bad to determine the dca
00549           if ((valid < (0.0004*rsq)) && (tries < 4))
00550              {batv.setX(pBach.x()*s2+xOut[i]);
00551               batv.setY(pBach.y()*s2+yOut[i]);
00552               batv.setZ(pBach.z()*s2+bachGeom2->origin().z());
00553               //Beginning of block 5.
00554               s1=(pBach.x()*dv0dotdb-dpV0.x())*diffc.x() + (pBach.y()*dv0dotdb-dpV0.y())*diffc.y() + (pBach.z()*dv0dotdb-dpV0.z())*diffc.z();
00555               s1=-s1/denom;
00556               v0atv.setX(dpV0.x()*s1+xV0.x());
00557               v0atv.setY(dpV0.y()*s1+xV0.y());
00558               v0atv.setZ(dpV0.z()*s1+xV0.z());
00559               //End of block 5 and beginning of block 5bis.
00560               //Cut: remove if V0 points away from Xi vertex
00561               if (((xV0.x()-v0atv.x())*pV0.x() + (xV0.y()-v0atv.y())*pV0.y() + (xV0.z()-v0atv.z())*pV0.z()) <= 0.) continue;
00562               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());
00563               //Cut: dca Xi daughters
00564               if (dca > (parsXi->dca_max*parsXi->dca_max)) continue;
00565               xpp.setX((v0atv.x()+batv.x())/2.);
00566               xpp.setY((v0atv.y()+batv.y())/2.);
00567               xpp.setZ((v0atv.z()+batv.z())/2.);
00568               //Cut: Xi decay length
00569               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;
00570               //Calculate xi impact parameter
00571               pXi.setX(pV0.x()+xOrig.x());
00572               pXi.setY(pV0.y()+xOrig.y());
00573               pXi.setZ(pV0.z()+xOrig.z());
00574               //Cut: remove if Xi points away from primary vertex
00575               if (((xpp.x()-xPvx.x())*pXi.x()+(xpp.y()-xPvx.y())*pXi.y()+(xpp.z()-xPvx.z())*pXi.z()) < 0.0) continue;
00576               //Calculate pt-Armenteros
00577               ptot_v02=ptV0*ptV0+pV0.z()*pV0.z();
00578               bdotx=xOrig.x()*pXi.x()+xOrig.y()*pXi.y()+xOrig.z()*pXi.z();
00579               vdotx=pV0.x()*pXi.x()+pV0.y()*pXi.y()+pV0.z()*pXi.z();
00580               if (bachGeom->charge() > 0)
00581                  {ppar=bdotx/pXi.mag();
00582                   pper=(ptot[k]*ptot[k]-ppar*ppar);
00583                   pper=(pper>0)? ::sqrt(ptot[k]*ptot[k]-ppar*ppar):0;
00584                  }
00585                   else
00586                  {ppar=vdotx/pXi.mag();
00587                   pper=(ptot_v02-ppar*ppar);
00588                   pper=(pper>0)? ::sqrt(ptot_v02-ppar*ppar):0;
00589                  }
00590               //Cut: pt-Armanteros
00591               if (pper > 0.33) continue;
00592               //Function helixDCA(charge,xpp,pXi,bxi);
00593               //helixDCA is defined in exi_c_utils.cc (pams/global/exi/).
00594               pt_tmp = ::sqrt(pXi.x()*pXi.x()+pXi.y()*pXi.y());
00595               bcharge_tmp = charge*(Bfield/kilogauss); //Field wanted in kGauss
00596               curvature_tmp = TMath::Abs(bcharge_tmp)*C_D_CURVATURE/pt_tmp;
00597               dip_tmp = atan(pXi.z()/pt_tmp);
00598               h_tmp = ((bcharge_tmp > 0) ? -1 : 1);
00599               phase_tmp = atan2(pXi.y(),pXi.x())-(h_tmp*C_PI_2);
00600               origin_tmp.setX(xpp.x());
00601               origin_tmp.setY(xpp.y());
00602               origin_tmp.setZ(xpp.z());
00603               StHelixD *globHelix = new StHelixD(curvature_tmp, dip_tmp, phase_tmp, origin_tmp, h_tmp);
00604               bxi=globHelix->distance(xPvx);
00605               //End of helixDCA
00606               //Cut: dca Xi to primary vertex
00607               if (bxi > parsXi->bxi_max)
00608                  {delete globHelix;
00609                   globHelix=0;
00610                   continue;
00611                   }
00612               StThreeVectorD p1 = globHelix->at(globHelix->pathLength(xPvx));
00613               StThreeVectorD p2(p1.x()-globHelix->xcenter(),p1.y()-globHelix->ycenter(),0);
00614               StThreeVectorD p3(xPvx.x()-globHelix->xcenter(),xPvx.y()-globHelix->ycenter(),0);
00615               if (p3.mag2() > p2.mag2()) bxi=-bxi;
00616               delete globHelix;
00617               globHelix=0;
00618               xiVertex = new StXiVertex();
00619               xiVertex->setPosition(xpp);
00620               xiVertex->addDaughter(trk[k]);
00621               xiVertex->setDcaBachelorToPrimaryVertex(trk[k]->impactParameter());
00622               xiVertex->setMomentumOfBachelor(xOrig);
00623               xiVertex->setDcaDaughters(::sqrt(dca));
00624               xiVertex->setDcaParentToPrimaryVertex(bxi);
00625               xiVertex->setV0Vertex(v0Vertex);
00626               
00628               //Set chi2 to trace SVT usage
00629               long int v0ChiSq = -1*(long int)v0Vertex->chiSquared();
00630               if(detId[k]==2 || detId[k]==3)
00631                  {//if an SVT track was used on bachelor set the last bit to 1
00632                   v0ChiSq |=((long int)1 << 0);
00633                   }
00634               v0ChiSq *=-1;
00635               xiVertex->setChiSquared((float)v0ChiSq);
00637               
00638               xiVertices.push_back(xiVertex);
00639               usedV0 = kTRUE;
00640               } //End if (valid && tries<4)
00641           } //End loop over 2 intersection points
00642       } // k-Loop
00643   delete bachGeom;
00644   bachGeom=0;
00645   delete bachGeom2;
00646   bachGeom2=0;
00647 
00648   if (alaCand)
00649      {alaCand = kFALSE;
00650       goto AssignSign;
00651       }
00652 
00653   return usedV0;
00654 }
00655 //_____________________________________________________________________________
00656 // $Id: StXiFinderMaker.cxx,v 1.23 2008/04/03 19:58:36 fisyak Exp $
00657 // $Log: StXiFinderMaker.cxx,v $
00658 // Revision 1.23  2008/04/03 19:58:36  fisyak
00659 // move parameters initialization from Init into InitRun
00660 //
00661 // Revision 1.22  2005/07/19 22:10:14  perev
00662 // STARFPE
00663 //
00664 // Revision 1.21  2004/04/15 20:15:53  jeromel
00665 // Forgot one of them
00666 //
00667 // Revision 1.20  2004/04/02 08:57:34  faivre
00668 // Use actual TPT flag rather than "not ITTF" for TPT tracks.
00669 //
00670 // Revision 1.19  2004/02/05 16:52:28  faivre
00671 // Add cut pt-Armanteros > 0.33 ; small cleanup.
00672 //
00673 // Revision 1.18  2004/02/04 14:53:36  faivre
00674 // Remove bug introduced in previous version (confusion dca/bxi when copy-pasting from helixDCA).
00675 //
00676 // Revision 1.17  2004/02/04 14:24:58  faivre
00677 // Add sign of dcaXiDaughters, slightly move cuts, small cleanup.
00678 //
00679 // Revision 1.16  2004/02/03 14:29:36  faivre
00680 // Spring-cleaning 2 months early ;-)  Algo strictly equivalent to previous one although huge reshaping.
00681 //
00682 // Revision 1.15  2004/02/02 12:10:54  faivre
00683 // XiFinder now able to run on muDsts as well :-)   + update user-friendliness.
00684 //
00685 // Revision 1.14  2003/11/08 18:26:00  faivre
00686 // Bfield + consistency int/short
00687 //
00688 // Revision 1.13  2003/09/17 12:00:23  faivre
00689 // RH8 : initialize everything in constructor.
00690 //
00691 // Revision 1.12  2003/09/02 17:58:59  perev
00692 // gcc 3.2 updates + WarnOff
00693 //
00694 // Revision 1.11  2003/08/22 17:47:14  caines
00695 // Get sign AND magnitude of mag field correctly for Xi and V0 finder
00696 //
00697 // Revision 1.10  2003/07/15 17:39:16  faivre
00698 // Doesn't take Bfield from gufld anymore (takes Bfield calculated in V0Finder).
00699 //
00700 // Revision 1.9  2003/07/04 17:52:46  faivre
00701 // Use SVT cuts if any dg has a SVT hit.
00702 //
00703 // Revision 1.8  2003/06/24 16:20:11  faivre
00704 // Uses SVT tracks. Fixed bool calculations. Exits when bad param. Reshaping.
00705 //
00706 // Revision 1.7  2003/05/14 19:15:36  faivre
00707 // Fancy choices Fortran/C++ V0's and Xi's. Xi rotating and like-sign.
00708 //
00709 // Revision 1.6  2003/05/07 10:51:42  faivre
00710 // Use brand new StHelixModel::setMomentum to solve memory leaks.
00711 //
00712 // Revision 1.4  2003/05/02 21:21:08  lbarnby
00713 // Now identify ITTF tracks by fittingMethod() equal to  kITKalmanFitId
00714 //
00715 // Revision 1.3  2003/04/30 20:38:22  perev
00716 // Warnings cleanup. Modified lines marked VP
00717 //
00718 // Revision 1.2  2003/04/30 19:16:04  faivre
00719 // Fix storage part. ITTF vs TPT Xis.
00720 //
00721 //