StRoot: StSecondaryVertexMaker/StKinkMaker.cxx Source File

00001 
00010 #include <iostream>
00011 #include <cstdlib>
00012 #include <cstring>
00013 #include "StMemStat.h"
00014 #include "StKinkMaker.h"
00015 #include "StKinkLocalTrack.hh"
00016 #include "StTrackGeometry.h"
00017 #include "StV0FinderMaker.h"
00018 #include "StEvent/StEventTypes.h"
00019 #include "StEvent.h"
00020 #include "TMath.h"
00021 #include "StTrack.h"
00022 #include "tables/St_tkf_tkfpar_Table.h"
00023 #include "StMessMgr.h"
00024 
00025 #include "math_constants.h"
00026 #include "phys_constants.h"
00027 #include "TVector2.h"
00028 #include "StThreeVectorF.hh"
00029 #include "TObjArray.h"
00030 #include "SystemOfUnits.h"
00031 
00032 
00033 #if !defined(ST_NO_NAMESPACES)
00034 using namespace units;
00035 #endif
00036 
00037 const double kaonMass= M_KAON_PLUS;
00038 const double pionMass= M_PION_PLUS;
00039 const double muonMass= M_MUON_PLUS;
00040 const double pi0Mass= M_PION_0;
00041 
00042 const double kaonToMuonQ= 0.236;
00043 const double kaonToPionQ= 0.205;
00044 const double pionToMuonQ= 0.030;
00045   
00046 const int MAXNUMOFTRACKS= 10000;
00047 const int SIZETRKIDCHECK= 1000;
00048 
00049 ClassImp(StKinkMaker)
00050 
00051 //=======================================================================
00052   StKinkMaker::StKinkMaker(const char *name):StMaker(name),m_tkfpar(0)
00053 {
00054   mTrack1     = 0;          
00055   mTrack2     = 0;   
00056   mGlobalTrks = 0;
00057   mParentTrackCandidate=0;             
00058   mDaughterTrackCandidate=0;
00059   mDaughterTrackUnic=0;
00060   mUseTracker = kTrackerUseBOTH;
00061   event       = 0;
00062   kinkVertex  = 0;
00063   mBfield     = -2.;   // arbitrary value, normalized after ...neat (gene)
00064 }
00065 //=======================================================================
00066 StKinkMaker::~StKinkMaker(){
00067 }
00068 
00069 //======================================================================
00070 
00079 Int_t StKinkMaker::Init(){
00080 
00081   // m_Mode -> SetTrackerUsage()
00082   if      (m_Mode == 1) SetTrackerUsage(kTrackerUseTPT);
00083   else if (m_Mode == 2) SetTrackerUsage(kTrackerUseITTF);
00084   else if (m_Mode == 3) SetTrackerUsage(kTrackerUseBOTH);
00085 
00086  
00087   return StMaker::Init();
00088 }
00089 //=============================================================================
00090 Int_t StKinkMaker::InitRun(int runumber) {
00091   m_tkfpar = (St_tkf_tkfpar*) GetDataBase("Calibrations/tracker/tkf_tkfpar");
00092   if (!m_tkfpar) {
00093     gMessMgr->Error(
00094       "StKinkMaker::InitRun() : could not find tkf_tkfpar in database.");
00095     return kStErr;
00096   }
00097   return StMaker::InitRun(runumber);
00098  }
00099 //=============================================================================
00100 Int_t StKinkMaker::Make(){//called for each event
00101 
00102 
00103   //******* variables
00104   unsigned short nNodes,i,j;
00105   int cutLast=0;
00106   StKinkLocalTrack* tempTrack;
00107   TObjArray trackArray(MAXNUMOFTRACKS);trackArray.SetOwner();
00108  
00109   tkf_tkfpar_st *tkfpar = m_tkfpar->GetTable();
00110   gMessMgr->Info()<<"StKinkMaker: impact parameter"<<tkfpar->impactCut<<endm;
00111   StThreeVectorF p,start,end;
00112 
00113   //******* GET 
00114   //### event
00115   event = (StEvent*)GetInputDS("StEvent");
00116   if (!event){
00117     gMessMgr->Warning("StKinkMaker:no StEvent;skip event");
00118     return kStWarn;
00119   }
00120 
00121   //### primary vertex
00122   StPrimaryVertex* vrtx=event->primaryVertex();
00123   if (!vrtx){
00124     gMessMgr->Warning("StKinkMaker:no primary vertex;skip event");
00125     return kStWarn;
00126   }
00127   mEventVertex = vrtx->position(); //StThreeVectorD
00128  
00129   //### global tracks to use
00130   StSPtrVecTrackNode& theNodes = event->trackNodes();
00131   nNodes = theNodes.size();
00132   mGlobalTrks=0;
00133   for (i=0;i<nNodes;i++){
00134     int nj = theNodes[i]->entries(global);
00135     for (j=0; j<nj;j++){
00136       StTrack* trk = theNodes[i]->track(global,j);
00137       if( !acceptTrack(trk) )continue;
00138 
00139      //******* find the magnetic field
00140       StTrackGeometry* trkGeom = trk->geometry();
00141       if(!trkGeom) continue; //ignore the track if it has no geometry 
00142 
00143       if(!mGlobalTrks)//calculate mBfield only for the first global trk from the event
00144         {
00145 
00146          StThreeVectorD p11 = trk->geometry()->momentum();
00147          StThreeVectorD p22 = trk->geometry()->helix().momentum(mBfield);
00148 
00149          if      (fabs(p22.x()) > 1.e-20) mBfield *= p11.x()/p22.x();
00150          else if (fabs(p22.y()) > 1.e-20) mBfield *= p11.y()/p22.y();
00151          else continue;
00152 
00153          if (fabs(mBfield) < 1.e-20) return kStWarn;
00154         }
00155         mGlobalTrks++;
00156       //### cut: fiducial volume
00157 
00158       Float_t trkStartRadius2D = trk->geometry()->origin().perp();
00159       Float_t trkEndRadius2D = trk->outerGeometry()->origin().perp();
00160       if (trkStartRadius2D < tkfpar->vertexRMin2D && //cut [133,179]
00161             trkEndRadius2D > tkfpar->vertexRMax2D ) 
00162         { continue;}
00163       if (trkStartRadius2D < tkfpar->vertexRMax2D ||
00164             trkEndRadius2D > tkfpar->vertexRMin2D ){
00165         tempTrack = new StKinkLocalTrack(trk);
00166         trackArray.Add(tempTrack);//keep the track if ok
00167 
00168       }
00169       }//for each track in the node
00170   }   //for each node
00171   //***************************************************************************************
00172   // cout<<" magnetic field="<<mBfield<<endl;
00173   //******* erase existing kinks; because everything runs in chain, the tkf makes its own kinks
00174   StSPtrVecKinkVertex& kinkVertices = event->kinkVertices();
00175   /* StMemStat memStat("kinkVertices");
00176   kinkVertices.getEntries();
00177   cout<<"memory used"<<  memStat.Used();
00178   */
00179   kinkVertices.erase(kinkVertices.begin(),kinkVertices.end());
00180   // kinkVertices.getEntries();
00181   //cout<<"memory used"<< memStat.Used();
00182   //  StSPtrVecKinkVertex kinkVertices2;
00183   //kinkVertices = kinkVertices2;
00184 
00185 //******* sorts by the trkStartRadius2D==>the potential parent trk, with 
00186       //smaller radius will be ahead of the potential daughters
00187   if (trackArray.GetEntries() > 0) trackArray.Sort();
00188   Int_t   kinkCandidate=0;
00189   Int_t cutPPt=0,cutPImpact=0,initial=0;
00190   
00191   int ni=trackArray.GetEntries();
00192   for (i=0;i<ni;i++){//parent
00193     initial++;
00194     mTrack1 = (StKinkLocalTrack*)trackArray.At(i);
00195     mParentTrackCandidate = mTrack1->trackBack();
00196     StTrackGeometry* myParentGeometry1 = mParentTrackCandidate->geometry(); //first point
00197     StTrackGeometry* myParentGeometry11 = mParentTrackCandidate->outerGeometry();//last point
00198     StPhysicalHelixD parentHelix = myParentGeometry1->helix();
00199     double parentPtot = myParentGeometry11->momentum().mag();
00200     double parentPt = myParentGeometry11->momentum().perp();
00201               
00202     //### cut pT parent minimum = 0.2 Gev/c
00203     if (myParentGeometry1->momentum().perp() < tkfpar->parentPtMin ) continue;
00204     cutPPt++;        
00205 
00206     //### cut : impact parameter parent
00207     mParentImpact = parentHelix.distance(mEventVertex);
00208     if (mParentImpact > tkfpar->impactCut ) continue;
00209     cutPImpact++;
00210     Int_t foundDaughters = 0;
00211     int nj = trackArray.GetLast()+1;
00212     for (j=i+1;j<nj;j++ ){//daughter
00213       mTrack2 = (StKinkLocalTrack*)trackArray.At(j);
00214       mDaughterTrackCandidate = mTrack2->trackBack();
00215       StTrackGeometry* myDaughterGeometry1 = mDaughterTrackCandidate->geometry();//first point
00216       StPhysicalHelixD daughterHelix = myDaughterGeometry1->helix();
00217       double daughterPtot = myDaughterGeometry1->momentum().mag();
00218       double daughterPt = myDaughterGeometry1->momentum().perp();
00219 
00220       //### cut: same charge
00221       if (myParentGeometry1->charge() != myDaughterGeometry1->charge() ) continue;
00222           
00223       //### cut:low lim fiducial volume = 133cm 
00224       if (myDaughterGeometry1->origin().perp() < tkfpar->vertexRMin2D ) continue; 
00225       
00226       //### cut:daughter impact parameter >2cm
00227       mDaughterImpact = daughterHelix.distance(mEventVertex);
00228       if (mDaughterImpact < mParentImpact) continue;
00229                 
00230       //### cut:last_pt-first_pt < 14 (radial) and < 20 (z)
00231       if (fabs(myDaughterGeometry1->origin().perp() - myParentGeometry11->origin().perp())
00232             > tkfpar->parentLastDaughterStart2D ) continue; 
00233       if (fabs(myDaughterGeometry1->origin().z() - myParentGeometry11->origin().z())
00234             > tkfpar->parentLastDaughterStartZ ) continue;
00235       cutLast++;
00236           
00237 
00238       //##############################################################
00239       //2D method (gets 3D aproximation after two 2D solution      //
00240       //    ***** START DETERMINATION 2D DCA .... good luck        //
00241      //###############################################################
00242       //kink vertex; look at the intersection pts in 2d of the 2helices
00243       pairD paths,path2,paths1;                             // helix pathLength vars;
00244       TVector2 r1,r2,xc1 ,xc2,tmp2V;                  // for 2D helices projections
00245       StThreeVectorD x1,x2,p1,p2,x1i,x2j ;//position and momentum at 2d dca
00246       double separation,dxc;    // helix circle params
00247       double dca_12,dca_12tmp; // 2D dca ... helpers
00248       double rad_1,rad_2;//radii of the projected helices
00249       double distanceOne,distanceTwo;
00250   
00251       //### parent
00252       xc1.Set(parentHelix.xcenter(),parentHelix.ycenter());//center of parent helix circle
00253       r1.Set(parentHelix.origin().x(),parentHelix.origin().y());//first pt of the helix
00254       r1 -= xc1;//distance
00255       rad_1 = r1.Mod();//sqrt(x*x+y*y) length of the radius
00256       //### daughter
00257       xc2.Set(daughterHelix.xcenter(),daughterHelix.ycenter());//center of daughter helix circle
00258       r2.Set(daughterHelix.origin().x(),daughterHelix.origin().y());//first pt on the helix
00259       r2 -= xc2;//distance
00260       rad_2 = r2.Mod();//sqrt(x*x+y*y) length of the radius
00261 
00262       tmp2V = xc1 - xc2;//distance between centers
00263       dxc = tmp2V.Mod();//length =O1O2
00264       separation = dxc - (rad_1+rad_2);//O1O2-(r1+r2)           
00265      
00266       distanceOne = 0;
00267       dca_12 = -9999;
00268       dca_12tmp = -9999;
00269       if (dxc==0)continue;  
00270       if (separation < 0){
00271         if (dxc <= TMath::Abs(rad_1-rad_2)){//one helix circle, completely inside the other
00272             double why = xc1.Y()-((rad_1+rad_2+dxc)*0.5 * (xc1.Y()-xc2.Y()))/dxc;
00273             double ecs = xc1.X()-((rad_1+rad_2+dxc)*0.5 * (xc1.X()-xc2.X()))/dxc;
00274             paths.first  = parentHelix.pathLength(ecs,why);
00275             paths.second = daughterHelix.pathLength(ecs,why);
00276             x1 = parentHelix.at(paths.first);
00277             x2 = daughterHelix.at(paths.second);
00278             dca_12=x1.z()-x2.z();
00279             }else {//2 intersection points=>2 solution, process those who aren't nan
00280                    //paths containes the path lengths for this solution with
00281                     //that for track 1 stores in 'first', and the track2 stored in 'second'
00282                   path2 = parentHelix.pathLength(rad_2,xc2.X(),xc2.Y());
00283                   if (!isnan(path2.first)){
00284                       paths.first = path2.first;
00285                       x1 = parentHelix.at(paths.first);
00286                       paths.second = daughterHelix.pathLength(x1.x(),x1.y());
00287                       x2 = daughterHelix.at(paths.second);
00288                       dca_12 = x1.z()-x2.z();
00289                       distanceOne = (x1-myDaughterGeometry1->origin()).mag();
00290                   }
00291                   if ((!isnan(path2.second))&&(path2.second!=path2.first))
00292                        { paths1.first = path2.second;
00293                          x1i = parentHelix.at(paths1.first);
00294                          paths1.second = daughterHelix.pathLength(x1i.x(),x1i.y());
00295                          x2j = daughterHelix.at(paths1.second);
00296                          distanceTwo = (x2j-myDaughterGeometry1->origin()).mag();
00297                          dca_12tmp = x1i.z()-x2j.z();
00298 
00299                          if( distanceTwo < distanceOne){
00300 
00301                       //second solution is better
00302                          dca_12 = dca_12tmp;
00303                          x1 = x1i;
00304                          x2 = x2j;
00305                          paths = paths1;
00306                       }//distance
00307                   
00308                    }else if (isnan(path2.second))//no solution
00309                       { continue;}
00310             }
00311             
00312        }else if(separation==0){
00313             double why = xc1.Y()-(rad_1 * (xc1.Y()-xc2.Y()))/dxc;
00314             double ecs = xc1.X()-(rad_1 * (xc1.X()-xc2.X()))/dxc;
00315             paths.first  = parentHelix.pathLength(ecs,why);
00316             paths.second = daughterHelix.pathLength(ecs,why);
00317             x1 = parentHelix.at(paths.first);
00318             x2 = daughterHelix.at(paths.second);
00319             dca_12=x1.z()-x2.z();
00320        }else if ((separation < tkfpar->dcaParentDaughterMax)){
00321               //helix circles are close, but not overlapping
00322               //find dca to point halfway between circle centers
00323                   tmp2V = (xc1 + xc2) * 0.5;
00324                   paths.first  = parentHelix.pathLength(tmp2V.X(),tmp2V.Y());
00325                   paths.second = daughterHelix.pathLength(tmp2V.X(),tmp2V.Y());
00326                   x1 = parentHelix.at(paths.first);
00327                   x2 = daughterHelix.at(paths.second);
00328                   dca_12=x1.z()-x2.z();
00329                   //   cout << "3rd case :xtarget = " << x1.x() << "\ty target =" << x1.y()<<endl;
00330                   }else {continue ;}//helix circle too far apart 
00331      
00332 //### cut: kink vertex in fiducial volume [133,179]
00333       if ((x1.perp() > tkfpar->vertexRMax2D)|| (x1.perp() < tkfpar->vertexRMin2D) )  continue;
00334       if ((x2.perp() > tkfpar->vertexRMax2D)|| (x2.perp() < tkfpar->vertexRMin2D) ) continue;
00335 //### cut: projectPointZDiff (2cm)
00336      if (fabs(dca_12) > tkfpar->projectPointZDiff) continue;
00337     
00338 //at this point, dca_12 can be used in 3D calculation of the dca
00339      //###############################################################
00340      //DONE @D DCA < START 3D DCA FINDING#                         //
00341      //##############################################################
00342 
00343           //### momentum of the 2 tracks , at the inters pt
00344      mParentMoment = parentHelix.momentumAt(paths.first, mBfield);
00345      mDaughterMoment = daughterHelix.momentumAt(paths.second, mBfield);
00346           
00347           //### cut: decay angle > 1
00348       mDecayAngle = (1./degree) * mParentMoment.angle(mDaughterMoment);
00349       if (mDecayAngle<tkfpar->thetaMin) continue;
00350   
00351 //### 3D dca between parent and daughter (comes out square)
00352      StThreeVectorD temp3V;
00353      double cos12,sin2_12,t1,t2;                     // 3D dca calculation vars
00354      double temp;
00355   
00356      p1 = mParentMoment/parentPtot;
00357      p2 = mDaughterMoment/daughterPtot;
00358 
00359      cos12 = p1.dot(p2);
00360      sin2_12 = (1.0 - cos12)*(1.- cos12);
00361      if( sin2_12){
00362         temp = dca_12/sin2_12;
00363         t1 = (-p1.z()+p2.z()*cos12)*temp;
00364         t2 = ( p2.z()-p1.z()*cos12)*temp;
00365 
00366         temp = rad_1*(parentPtot/parentPt);
00367         temp *= sin(t1/temp);
00368         x1i = x1 + p1.pseudoProduct(temp,temp,t1);
00369 
00370         temp = rad_2*(daughterPtot/daughterPt);
00371         temp *= sin(t2/temp);
00372         x2j = x2 + p2.pseudoProduct(temp,temp,t2);
00373 
00374         dca_12tmp = (x1i - x2j).mag2();
00375         dca_12 *= dca_12;
00376 
00377         if( dca_12tmp < dca_12){
00378             paths.first  = parentHelix.pathLength(x1i.x(),x1i.y());
00379             paths.second = daughterHelix.pathLength(x2j.x(),x2j.y());
00380             x1i =  parentHelix.at(paths.first);
00381             x2j =  daughterHelix.at(paths.second);
00382               dca_12tmp = (x1i - x2j).mag2();
00383             if( dca_12tmp < dca_12){
00384                 x1 = x1i;
00385                 x2 = x2j;
00386                 mParentMoment = parentHelix.momentumAt(paths.first, mBfield);
00387                 mDaughterMoment = daughterHelix.momentumAt(paths.second,mBfield);
00388                 dca_12 = dca_12tmp;
00389                 }
00390         }
00391      }
00392      
00393 
00394   //############ end 3D DCA CALCUALTION
00395    mDca = sqrt(dca_12);
00396 
00397 
00398   //###cut: 3D dca < .5
00399   if (mDca>(tkfpar->dcaParentDaughterMax) ) continue;
00400   // cout << " dca is = " << mDca;       
00401   //cutDca++;
00402 
00403   mKinkVertex = (x1 + x2) * 0.5;//vertex position
00404   Float_t dx,dy;
00405   dx = mKinkVertex.x() - myParentGeometry11->origin().x();
00406   dy = mKinkVertex.y() - myParentGeometry11->origin().y();
00407   Float_t distanceKinkParent2D   = sqrt(dx*dx+dy*dy );
00408   dx = mKinkVertex.x() - myDaughterGeometry1->origin().x();
00409   dy = mKinkVertex.y() - myDaughterGeometry1->origin().y();
00410   Float_t distanceKinkDaughter2D = sqrt(dx*dx+dy*dy );
00411   Float_t distanceKinkParentZ    = mKinkVertex.z() - myParentGeometry11->origin().z();
00412   Float_t distanceKinkDaughterZ  = mKinkVertex.z() - myDaughterGeometry1->origin().z();
00413           
00414   //cut: distance kink track :radial<14, z<20
00415   if (distanceKinkParent2D > tkfpar->distanceKinkParent2D ) continue; 
00416   if (distanceKinkDaughter2D > tkfpar->distanceKinkDaughter2D ) continue; 
00417   if (distanceKinkParentZ > tkfpar->distanceKinkParentZ ) continue; 
00418   if (distanceKinkDaughterZ > tkfpar->distanceKinkDaughterZ ) continue; 
00419   foundDaughters++;
00420   mDaughterTrackUnic = myDaughterGeometry1;
00421    }//loop j.. daughters
00422   if(foundDaughters!=1)continue;//if found more than 1 daugthers for one parent track,ignore it
00423 
00424   kinkCandidate++;                
00425   FillEvent(mDaughterTrackUnic,myParentGeometry11);
00426   kinkVertices.push_back(kinkVertex);
00427   //  cout<<"foundDaugthers"  <<foundDaughters<<endl;
00428 
00429  }//loop i .. parents
00430   // trackArray.Delete(); 
00431   gMessMgr->Info() << "StKinkMaker:: Found " << kinkCandidate << " kink candidates " << endm;
00432 
00433   //========================================================================
00434   // Look for kinks in which 2 different parents are sharing the same daughter
00435   // mark the coresonding  kink vertices by changing the decayAngle value
00436   // the new decay angle is [old*100+99]
00437   // make the vertex found zombie ???!!!
00438 
00439   if(kinkVertices.size()>1) Crop();
00440   return kStOK; 
00441 }
00442         
00443 //######################### DONE!!! ###############################################
00444      
00445 
00447 void StKinkMaker::FillEvent(StTrackGeometry *myDaughterGeometry1,StTrackGeometry *myParentGeometry11){
00448 
00449   kinkVertex = new StKinkVertex();
00450   StThreeVectorF pMomMinusdMom = mParentMoment - mDaughterMoment;
00451   Float_t  deltaKaonMuon = fabs(sqrt(mParentMoment.mag2() + kaonMass*kaonMass)   -
00452                                 sqrt(mDaughterMoment.mag2() + muonMass*muonMass) - 
00453                                 pMomMinusdMom.mag());
00454   Float_t  deltaKaonPion = fabs(sqrt(mParentMoment.mag2() + kaonMass*kaonMass)   -
00455                                 sqrt(mDaughterMoment.mag2() + pionMass*pionMass) - 
00456                                 sqrt(pMomMinusdMom.mag2() + pi0Mass*pi0Mass));
00457   Float_t  deltaPionMuon = fabs(sqrt(mParentMoment.mag2() + pionMass*pionMass)   -
00458                                 sqrt(mDaughterMoment.mag2() + muonMass*muonMass) - 
00459                                 pMomMinusdMom.mag());  
00460   //K-/+=>pi+/- + pi0 
00461   if ((deltaKaonPion < deltaKaonMuon) && (deltaKaonPion < deltaPionMuon) ){
00462       Float_t asinArg = (mDaughterMoment.mag() / kaonToPionQ) * sin(mDecayAngle*degree);//sin(theta_cm)
00463       if (fabs(asinArg) < 1. ) {
00464           kinkVertex->setDecayAngleCM(float((1./degree) * asin(asinArg)));
00465         } 
00466       else {
00467           kinkVertex->setDecayAngleCM(999.) ;
00468         }
00469       if( myParentGeometry11->charge() > 0 ){ //parent K+
00470           kinkVertex->setGeantIdDaughter(8);//daughter is pi+
00471           kinkVertex->setGeantIdParent(11);
00472         }
00473       else {// parent K-
00474           kinkVertex->setGeantIdDaughter(9);//daughter pi-
00475           kinkVertex->setGeantIdParent(12);
00476         }
00477     } //K-/+=>mu+/- + (anti)neutrino
00478   else if ((deltaKaonMuon < deltaKaonPion) && (deltaKaonMuon < deltaPionMuon) ) {
00479       Float_t asinArg = (mDaughterMoment.mag() / kaonToMuonQ) * sin(mDecayAngle*degree);
00480       if (fabs(asinArg) < 1. ) {
00481           kinkVertex->setDecayAngleCM( (1./degree) * asin(asinArg) );
00482         } 
00483       else {
00484           kinkVertex->setDecayAngleCM( 999.) ;
00485         }
00486       if( myParentGeometry11->charge() > 0 ){ //parent K+
00487           kinkVertex->setGeantIdDaughter(5);//daughter mu +
00488           kinkVertex->setGeantIdParent(11);
00489         } 
00490       else {//parent K-
00491           kinkVertex->setGeantIdDaughter(6);//daughter mu-
00492           kinkVertex->setGeantIdParent(12);
00493         }   
00494     } else { //pi+/-=>mu+/- + (anti)neutrino
00495         Float_t asinArg = (mDaughterMoment.mag() / pionToMuonQ) * sin(mDecayAngle*degree);
00496         if( fabs(asinArg) < 1. ) {
00497             kinkVertex->setDecayAngleCM( (1./degree) * asin(asinArg) );  
00498           } 
00499         else {
00500             kinkVertex->setDecayAngleCM( 999.) ;
00501           }
00502         if (myParentGeometry11->charge() > 0 ){ //parent pi+
00503             kinkVertex->setGeantIdDaughter(5);//daughter mu+
00504             kinkVertex->setGeantIdParent(8);
00505           } 
00506         else {//parent pi-
00507             kinkVertex->setGeantIdDaughter(6);//daughter mu-
00508             kinkVertex->setGeantIdParent(9);
00509           }   
00510       }
00511  
00512   kinkVertex->setDcaParentDaughter(mDca);
00513   kinkVertex->setDcaDaughterPrimaryVertex(mDaughterImpact);//dca daughter from event vertex
00514   kinkVertex->setDcaParentPrimaryVertex(mParentImpact);//dca parent from event vertex
00515   
00516   //distance last point parent - first pt daughter
00517   kinkVertex->setHitDistanceParentDaughter( (myDaughterGeometry1->origin() - myParentGeometry11->origin()).mag() );
00518         
00519   //distance last point parent primary vertex
00520   kinkVertex->setHitDistanceParentVertex(sqrt( pow(mKinkVertex.x() - myParentGeometry11->origin().x(), 2) +
00521                           pow(mKinkVertex.y() - myParentGeometry11->origin().y(), 2) +
00522                           pow(mKinkVertex.z() - myParentGeometry11->origin().z(), 2)) );
00523 
00524   //set dE-delta enrgy for different deay hypotheses
00525   kinkVertex->setdE(1,deltaKaonMuon);//dE for K->mu + (anti)neutrino
00526   kinkVertex->setdE(2,deltaKaonPion);//dE for K->pi+pi0
00527   kinkVertex->setdE(3,deltaPionMuon);//dE for pi->mu +(anti)neutrino
00528   
00529   //set momentum
00530   kinkVertex->setParentMomentum(mParentMoment);
00531   kinkVertex->setDaughterMomentum(mDaughterMoment);
00532   kinkVertex->setParent(mParentTrackCandidate);
00533 
00534   //set dcay angle        
00535   kinkVertex->setDecayAngle(mDecayAngle);
00536 
00537   //add daughter
00538   kinkVertex->addDaughter(mDaughterTrackCandidate);
00539  
00540   //set kinkVertex position
00541   kinkVertex->setPosition(mKinkVertex); 
00542 }
00543 
00544 
00559 bool StKinkMaker::acceptTrack(StTrack *trk)
00560 {
00561   //cout << "DEBUG Track [" << trk->fittingMethod() << "] [" << trk->flag() << "] FitIds [" 
00562   //     << kKalmanFitId << "] [" << kITKalmanFitId << "] selector [" << mUseTracker << "]" << endl;
00563 
00564   // cut on flag
00565   if (trk->flag() <= 0) return false;
00566 
00567   // on fittingMethod() 
00568   if  (  
00569        //( trk->fittingMethod() == kHelix3DIdentifier  && (mUseTracker == kTrackerUseTPT  || mUseTracker == kTrackerUseBOTH) ) ||
00570          ( trk->fittingMethod() != kITKalmanFitId  && (mUseTracker == kTrackerUseTPT  || mUseTracker == kTrackerUseBOTH) ) ||
00571          ( trk->fittingMethod() == kITKalmanFitId  && (mUseTracker == kTrackerUseITTF || mUseTracker == kTrackerUseBOTH) ) ){
00572     return true;
00573   } else {
00574     return false;
00575   }
00576 
00577 }
00578 
00579 
00596 void StKinkMaker::SetTrackerUsage(Int_t opt)
00597 {
00598   mUseTracker=opt;
00599   gMessMgr->Info() << "StKinkMaker::SetTrackerUsage : Setting option to " << mUseTracker << endm;
00600 }
00601                 
00602 void StKinkMaker::Crop(){
00603   // Loop over kinks and remove (makeZombie() )+ change the decay angle for those that have one daughter sharing two different parents
00604   gMessMgr->Info()<<"StKinkMaker::Crop() : Starting ...."<<endm;
00605   event = (StEvent *)GetInputDS("StEvent");
00606   StSPtrVecKinkVertex& kinkVertices = event->kinkVertices();
00607   StKinkVertex* kinkVtxPtr1 = 0;//new StKinkVertex();
00608   StKinkVertex* kinkVtxPtr2 = 0;// new StKinkVertex();
00609   StTrack* daughterTrk1 = 0;
00610   StTrack* daughterTrk2 = 0;
00611   int iKinks = kinkVertices.size();
00613   // take pairs of consecutive kink vertices and check for same track key() for daugther
00614   for(Int_t ikv1=0; ikv1<iKinks-1; ikv1++) {//first kink vertex from teh container
00615     kinkVtxPtr1 = kinkVertices[ikv1];
00616     daughterTrk1=kinkVtxPtr1->daughter();
00617     //  gMessMgr->Info()<<"###########StKinkMaker:daughter1->key()"<<daughterTrk1->key()<<endm;                  
00618     for(Int_t ikv2=ikv1+1; ikv2<iKinks; ikv2++) {//next kink vertex in the container
00619 
00620       kinkVtxPtr2 = kinkVertices[ikv2]; 
00621       daughterTrk2=kinkVtxPtr2->daughter();
00622     //  gMessMgr->Info()<<"###########StKinkMaker: daughter2->key()"<<daughterTrk2->key()<<endm;  
00623       if(daughterTrk1->key() == daughterTrk2->key()) {
00624         float decAngKinkVtx1 = kinkVtxPtr1->decayAngle();
00625         float decAngKinkVtx2 = kinkVtxPtr2->decayAngle(); 
00626         kinkVtxPtr1->setDecayAngle(decAngKinkVtx1*100+99);
00627         kinkVtxPtr2->setDecayAngle(decAngKinkVtx2*100+99);        
00628       }
00629 
00630     }//second vertex loop
00631   }//first vtx loop
00632 
00634   // make Zombie the kink vertices with problems
00635  //  for(int ikv=iKinks-1;ikv>=0;ikv--){
00636 //     kinkVertex = kinkVertices[ikv];
00637 //     if( (kinkVertex->geantIdDaughter()%100) == 99){
00638 
00639 //       kinkVertex->makeZombie();
00640 //       //   gMessMgr->Info()<<"###########StKinkMaker: AFTER "<<kinkVertex->decayAngle()%100<<endm;  
00641 //       iKinks--;
00642 //     }
00643 //   }
00644 
00645 //   gMessMgr->Info()<<"StKinkMaker::Crop() : saving "<< iKinks << " Kink Candidates" <<endm;
00646 
00647 
00648 }