/Home/aooliver/hepsoft/lcdroot/LCDRootApps/PhUtil/src/LCDVToplBase.cxx

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------------
00002 // $Id: LCDVToplBase.cxx,v 1.2 2001/10/08 17:53:59 toshi Exp $
00003 // ----------------------------------------------------------------------------
00004 
00005 #include "LCDVToplBase.h"
00006 
00007 #include <iostream.h>
00008 
00009 ClassImp(LCDVToplBase)
00010 
00011 //__________________________________________________________________________
00012 //Constructor and Deconstructor
00013 LCDVToplBase::LCDVToplBase() {
00014   InitBase();
00015   CalcIPErrorMatrix();
00016 }
00017 //__________________________________________________________________________
00018 
00019 //__________________________________________________________________________
00020 LCDVToplBase::~LCDVToplBase(){
00021   CleanBase(1);
00022 }
00023 //__________________________________________________________________________
00024 
00025 //__________________________________________________________________________
00026 void LCDVToplBase::SetUp(const TVector3& ipv, 
00027                          const TVector3& bfldv) {
00028 
00029   if (m_dblevel > 9) 
00030     cout << " <<LCDVToplBase::SetUp1>> Start..." << endl;
00031 
00032   //This loop to make sure position and momentum in trk at 2d POCA
00033   // and  to make LCDVToplTRK.
00034   TVector3 pja(0.0,0.0,0.0); 
00035   if (m_pja.Mag2() == 0 && m_nvtrack ) {
00036     for (Int_t i=0; i < m_nvtrack ; i++) {
00037       //Accumurate track momentum.
00038       pja += GetTrackAt(i)->GetMomentumVector(0.0);
00039     } //loop over zxtrks
00040   }
00041 
00042   SetUp(ipv,bfldv,pja);
00043 
00044   if (m_dblevel > 9) 
00045     cout << " <<LCDVToplBase::SetUp1>> End..." << endl;
00046 }
00047 //__________________________________________________________________________
00048 
00049 //__________________________________________________________________________
00050 void LCDVToplBase::SetUp(const TVector3& ipv, 
00051                          const TVector3& bfldv,
00052                          const TVector3& jeta) {
00053   if (m_dblevel > 9) 
00054     cout << " <<LCDVToplBase::SetUp2>> Start..." << endl;
00055 
00056   CleanBase(0);
00057 
00058   SetIP(ipv);
00059   SetBField(bfldv);
00060   m_pja=jeta;
00061   CalcIPErrorMatrix();
00062 
00063   if (m_dblevel > 9) 
00064     cout << " <<LCDVToplBase::SetUp2>> End..." << endl;
00065 }
00066 //__________________________________________________________________________
00067 
00068 //__________________________________________________________________________
00069 //Init...
00070 void LCDVToplBase::InitBase() {
00071   m_nvtrack       = 0;
00072   m_vertex_list   = new TClonesArray("LCDVToplVRT",100);
00073   m_nvertex       = 0;
00074   m_vtrack_list   = new TClonesArray("LCDVToplTRK",100);
00075   m_nvrttrack     = 0;
00076   for(Int_t i=0 ; i < 100 ; i++) {
00077     m_vrttrack_list[i]=0;
00078   }
00079     
00080   m_ripe  = 3.0e-4;
00081   m_zipe  = 6.0e-4;
00082 
00083   // Initialize jet total momentum
00084   m_pja.SetXYZ(0.0,0.0,0.0); 
00085 
00086   //Set Beam Position (Temporary)  
00087   m_ipv.SetXYZ(0.0,0.0,0.0); 
00088 
00089   // m_dblevel=0;
00090 
00091   //
00092   // Pt corrected mass relates.
00093   //
00094   m_Rmass    = 0.;
00095   m_PtMass   = 0.;
00096   m_PtMassEr = 0.;
00097   m_Masstag  = 0.;
00098   m_MissPt   = 0.;
00099 
00100   m_SeedVrtRMax  =2.3  ;
00101   m_SeedVrtK0Mass=0.015;
00102   m_SeedVrtLambda0Mass=0.015;
00103   m_SeedVrtSigMin=2.0  ;
00104 
00105   m_TrdiMax=0.1;
00106   m_LodiMin=0.05;
00107   m_LodrMin=0.25;
00108   m_LodrMax=100.0;
00109   m_AntaMax=TMath::Pi()/2.0;
00110 
00111   m_PTIP[0]=1.5*m_ripe;
00112   m_PTIP[1]=1.5*m_zipe;
00113   m_NVsigma=1.0;
00114 
00115   m_UseNN     =kFALSE;
00116   m_SeedVrtNN =0.4;
00117   m_TrackAttNN=0.6;
00118   m_SelectNN  =0.5;
00119 
00120   m_p4Vertex.SetXYZT(0.,0.,0.,0.);
00121 
00122   m_mcjetflav=0;
00123 
00124   m_ntrk2vrt=0;
00125   m_trk2vrt=0;
00126   m_nvrt2trk=0;
00127   m_vrt2trk=0;
00128   m_trkl=0;
00129 
00130 }
00131 //__________________________________________________________________________
00132 
00133 //__________________________________________________________________________
00134 void LCDVToplBase::DeleteTracks(Int_t iflg){
00135   if (m_vtrack_list)   m_vtrack_list->Delete();
00136   m_nvtrack=0;
00137   if (iflg) {
00138     delete m_vtrack_list;
00139     m_vtrack_list=0;
00140   }
00141   m_nvrttrack=0;
00142   for(Int_t i=0 ; i < 100 ; i++) {
00143     m_vrttrack_list[i]=0;
00144   }
00145 }
00146 //__________________________________________________________________________
00147 
00148 //__________________________________________________________________________
00149 void LCDVToplBase::DeleteVertices(Int_t iflg){
00150   if (m_vertex_list) m_vertex_list->Delete();
00151   if (iflg) {
00152     delete m_vertex_list;
00153     m_vertex_list=0;
00154   }
00155 }
00156 //__________________________________________________________________________
00157 
00158 //__________________________________________________________________________
00159 //Delete LCDVToplVRT, LCDVToplTRK...
00160 void LCDVToplBase::CleanBase(Int_t iflg) {
00161   DeleteTracks(iflg);
00162   DeleteVertices(iflg);
00163   if (m_trk2vrt) {
00164     delete [] m_trk2vrt;
00165     m_trk2vrt=0;
00166   }
00167   if (m_vrt2trk) {
00168     delete [] m_vrt2trk;
00169     m_vrt2trk=0;
00170   }
00171   if (m_trkl) {
00172     delete [] m_trkl;
00173     m_trkl=0;
00174   }
00175   m_ntrk2vrt=0;
00176   m_nvrt2trk=0;
00177 }
00178 //__________________________________________________________________________
00179 
00180 //__________________________________________________________________________
00181 void LCDVToplBase::CalcIPErrorMatrix(){
00182   Double_t ripe2=m_ripe*m_ripe;
00183   Double_t zipe2=m_zipe*m_zipe;
00184   m_iper[0]=ripe2/2.0;
00185   m_iper[1]=0.0      ; m_iper[2]=ripe2/2.0;
00186   m_iper[3]=0.0      ; m_iper[4]=0.0      ; m_iper[5]=zipe2;
00187   TMatrixD t(3,3);
00188   GetBeamIPErrorMatrix(t);
00189   Double_t det_iper;
00190   t.Invert(&det_iper);
00191   SetBeamIPErrorMatrixInv(t);
00192 }
00193 //__________________________________________________________________________
00194 
00195 //__________________________________________________________________________
00196 void LCDVToplBase::GetBeamIPErrorMatrix(TMatrixD& t) {
00197   t(0,0)=m_iper[0]; t(0,1)=m_iper[1]; t(0,2)=m_iper[3];
00198   t(1,0)=m_iper[1]; t(1,1)=m_iper[2]; t(1,2)=m_iper[4];
00199   t(2,0)=m_iper[3]; t(2,1)=m_iper[4]; t(2,2)=m_iper[5];
00200 }
00201 //__________________________________________________________________________
00202 
00203 //__________________________________________________________________________
00204 void LCDVToplBase::GetBeamIPErrorMatrixInv(TMatrixD& t) {
00205   t(0,0)=m_iper_inv[0]; t(0,1)=m_iper_inv[1]; t(0,2)=m_iper_inv[3];
00206   t(1,0)=m_iper_inv[1]; t(1,1)=m_iper_inv[2]; t(1,2)=m_iper_inv[4];
00207   t(2,0)=m_iper_inv[3]; t(2,1)=m_iper_inv[4]; t(2,2)=m_iper_inv[5];
00208 }
00209 //__________________________________________________________________________
00210 
00211 //__________________________________________________________________________
00212 void LCDVToplBase::SetBeamIPErrorMatrix(TMatrixD& t) {
00213   m_iper[0]=t(0,0);
00214   m_iper[1]=t(1,0); m_iper[2]=t(1,1);
00215   m_iper[3]=t(2,0); m_iper[4]=t(2,1); m_iper[5]=t(2,2);
00216 }
00217 //__________________________________________________________________________
00218 
00219 //__________________________________________________________________________
00220 void LCDVToplBase::SetBeamIPErrorMatrixInv(TMatrixD& t) {
00221   m_iper_inv[0]=t(0,0);
00222   m_iper_inv[1]=t(1,0); m_iper_inv[2]=t(1,1);
00223   m_iper_inv[3]=t(2,0); m_iper_inv[4]=t(2,1); m_iper_inv[5]=t(2,2);
00224 }
00225 //__________________________________________________________________________
00226 
00227 //__________________________________________________________________________
00228 LCDVToplTRK* LCDVToplBase::GetTrackAt(Int_t itrk) {
00229   return (LCDVToplTRK*)(GetTracks()->At(itrk));
00230 }
00231 //__________________________________________________________________________
00232 
00233 //__________________________________________________________________________
00234 Int_t LCDVToplBase::GetTrackIdx(LCDVToplTRK* trk) {
00235   return GetTracks()->IndexOf(trk);
00236 }
00237 //__________________________________________________________________________
00238 
00239 //__________________________________________________________________________
00240 //Get # of tracks
00241 Int_t LCDVToplBase::GetTrackEntries(){
00242   return m_nvtrack;
00243   // return GetTracks()->GetEntries();
00244 }
00245 //__________________________________________________________________________
00246 
00247 //__________________________________________________________________________
00248 LCDVToplVRT* LCDVToplBase::GetVertexAt(Int_t ivrt) {
00249   return (LCDVToplVRT*)(GetVertices()->At(ivrt));
00250 }
00251 //__________________________________________________________________________
00252 
00253 //__________________________________________________________________________
00254 Int_t LCDVToplBase::GetVertexIdx(LCDVToplVRT* vrt) {
00255   return GetVertices()->IndexOf(vrt);
00256 }
00257 //__________________________________________________________________________
00258 
00259 //__________________________________________________________________________
00260 Int_t LCDVToplBase::GetVertexEntries(){
00261   return m_nvertex;
00262   //    return GetVertices()->GetEntries();
00263 }
00264 //__________________________________________________________________________
00265 
00266 //__________________________________________________________________________
00267 Int_t LCDVToplBase::GetVertexEntriesExceptV0(){
00268   Int_t nvrt=GetVertexEntries();
00269   Int_t nvrtnov0=0;
00270   LCDVToplVRT* vrt;
00271   for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
00272     vrt=GetVertexAt(ivrt);
00273     if (CheckKs0Vertex(vrt))     continue;
00274     if (CheckLambda0Vertex(vrt)) continue;
00275     nvrtnov0++;
00276   }
00277   return nvrtnov0;
00278 }
00279 //__________________________________________________________________________
00280 
00281 //__________________________________________________________________________
00282 LCDVToplTRK* LCDVToplBase::AddTrack(Int_t tkid, 
00283                                     Double_t* tkpar, 
00284                                     Double_t* e_tkpar, 
00285                                     Double_t chrg, 
00286                                     Double_t mfld, 
00287                                     Int_t mc_ind, 
00288                                     const TVector3& ipv) {
00289   LCDVToplTRK* trk = new((*m_vtrack_list)[m_nvtrack++]) 
00290     LCDVToplTRK(tkid,tkpar,e_tkpar,chrg,mfld,mc_ind,ipv);
00291   m_pja += trk->GetMomentumVector(0.0);
00292   return trk;
00293 }
00294 //__________________________________________________________________________
00295 
00296 //__________________________________________________________________________
00297 //Get Initial vertexposition rrack and track 
00298 TVector3 LCDVToplBase::GetInitialVertexT2T(LCDVToplTRK* trk1, 
00299                                            LCDVToplTRK* trk2,
00300                                            Double_t l1, Double_t l2) {
00301   //initial vertex point for given tracks
00302 
00303   TVector3 posv1   = trk1->GetPositionVector(l1);
00304   TVector3 momv1   = trk1->GetMomentumVector(l1);
00305   TVector3 ip2pos1(m_ipv - posv1);
00306   if (l1 != 0.0) {
00307     Double_t mu1   = momv1.X()*ip2pos1.X() + momv1.Y()*ip2pos1.Y();
00308     mu1 /= (momv1.X()*momv1.X() + momv1.Y()*momv1.Y());
00309     posv1 += mu1*momv1;
00310   }
00311   // Double_t omega1  = trk1->GetTrackParameter(2);
00312   // Double_t phi0pr1 = trk1->GetTrackParameter(1) - TMath::Pi()/2.0 + l1/omega1;
00313   Double_t phi0pr1 = momv1.Phi() - TMath::Pi()/2.0;
00314   Double_t cspr1   = TMath::Cos(phi0pr1);
00315   Double_t snpr1   = TMath::Sin(phi0pr1);
00316   Double_t tanlpr1 = trk1->GetTrackParameter(4);
00317   Double_t z0pr1   = posv1.Z();
00318   Double_t x0pr1   = cspr1*(posv1[0] - m_ipv[0]) + snpr1*(posv1[1] - m_ipv[1]);
00319   Double_t s1spr1  = trk1->GetErrorTsiApprox(l1);
00320   Double_t s2spr1  = trk1->GetErrorEtaApprox(l1)*(1.0 + tanlpr1*tanlpr1);
00321 
00322   TVector3 posv2   = trk2->GetPositionVector(l2);
00323   TVector3 momv2   = trk2->GetMomentumVector(l2);
00324   TVector3 ip2pos2(m_ipv - posv2);
00325   if (l2 != 0.0) {
00326     Double_t mu2   = momv2.X()*ip2pos2.X() + momv2.Y()*ip2pos2.Y();
00327     mu2 /= (momv2.X()*momv2.X() + momv2.Y()*momv2.Y());
00328     posv2 += mu2*momv2;
00329   }
00330   // Double_t omega2  = trk2->GetTrackParameter(2);
00331   // Double_t phi0pr2 = trk2->GetTrackParameter(1) - TMath::Pi()/2.0 + l2/omega2;
00332   Double_t phi0pr2 = momv2.Phi() - TMath::Pi()/2.0;
00333   Double_t cspr2   = TMath::Cos(phi0pr2);
00334   Double_t snpr2   = TMath::Sin(phi0pr2);
00335   Double_t tanlpr2 = trk2->GetTrackParameter(4);
00336   Double_t z0pr2   = posv2.Z();
00337   Double_t x0pr2   = cspr2*(posv2[0] - m_ipv[0]) + snpr2*(posv2[1] - m_ipv[1]);
00338   Double_t s1spr2  = trk2->GetErrorTsiApprox(l2);
00339   Double_t s2spr2  = trk2->GetErrorEtaApprox(l2)*(1.0 + tanlpr2*tanlpr2);
00340 
00341   Double_t c1 = 
00342     cspr1*cspr1/s1spr1 + tanlpr1*tanlpr1*snpr1*snpr1/s2spr1 +
00343     cspr2*cspr2/s1spr2 + tanlpr2*tanlpr2*snpr2*snpr2/s2spr2;
00344   Double_t c2 =
00345     snpr1*snpr1/s1spr1 + tanlpr1*tanlpr1*cspr1*cspr1/s2spr1 +
00346     snpr2*snpr2/s1spr2 + tanlpr2*tanlpr2*cspr2*cspr2/s2spr2;
00347   Double_t c3 = 1/s2spr1 + 1/s2spr2;
00348   Double_t c4 =
00349     cspr1*snpr1/s1spr1 - tanlpr1*tanlpr1*snpr1*cspr1/s2spr1 +
00350     cspr2*snpr2/s1spr2 - tanlpr2*tanlpr2*snpr2*cspr2/s2spr2;
00351   Double_t c5 = tanlpr1*snpr1/s2spr1 + tanlpr2*snpr2/s2spr2;
00352   Double_t c6 =-(tanlpr1*cspr1/s2spr1 + tanlpr2*cspr2/s2spr2);
00353   Double_t c7 =
00354     cspr1*x0pr1/s1spr1 + z0pr1*tanlpr1*snpr1/s2spr1 +
00355     cspr2*x0pr2/s1spr2 + z0pr2*tanlpr2*snpr2/s2spr2;
00356   Double_t c8 =
00357     snpr1*x0pr1/s1spr1 - z0pr1*tanlpr1*cspr1/s2spr1 +
00358     snpr2*x0pr2/s1spr2 - z0pr2*tanlpr2*cspr2/s2spr2;
00359   Double_t c9 =z0pr1/s2spr1 + z0pr2/s2spr2;
00360 
00361   TMatrixD rot(3,3);
00362   rot(0,0) = c1; rot(0,1) = c4; rot(0,2)= c5;
00363   rot(1,0) = c4; rot(1,1) = c2; rot(1,2)= c6;
00364   rot(2,0) = c5; rot(2,1) = c6; rot(2,2)= c3;
00365   TMatrixD tmp(3,1);
00366   tmp(0,0) = c7;
00367   tmp(1,0) = c8;
00368   tmp(2,0) = c9;
00369   
00370   TVector3 xdv=m_ipv;
00371   Double_t det_rot;
00372   rot.Invert(&det_rot);
00373   if (det_rot) {
00374     TMatrixD tmp1(rot,TMatrixD::kMult,tmp);
00375     // tmp1.Mult(rot,tmp);
00376     xdv.SetXYZ(tmp1(0,0)+m_ipv.X(),tmp1(1,0)+m_ipv.Y(),tmp1(2,0));
00377     // xdv.SetXYZ(tmp1(0,0)+m_ipv[0],tmp1(1,0)+m_ipv[1],tmp1(2,0));
00378     // avoid being too far off
00379     // TVector3 xdv2(xdv[0],xdv[1],xdv[2]-m_ipv[2]);
00380     // if ( xdv2.Mag2() > 16.0 ) {
00381     //  xdv=m_ipv;
00382     // }
00383   }
00384 
00385   return xdv;
00386 }
00387 //__________________________________________________________________________
00388 
00389 //__________________________________________________________________________
00390 void LCDVToplBase::RemoveVertexAt(Int_t i) {
00391   RemoveVertex(GetVertexAt(i));
00392 }
00393 //__________________________________________________________________________
00394 
00395 //__________________________________________________________________________
00396 void LCDVToplBase::RemoveVertex(LCDVToplVRT* vrt) {
00397   GetVertices()->Remove(vrt);
00398   CompressVerticesList();
00399   m_nvertex--;
00400 }
00401 //__________________________________________________________________________
00402 
00403 //__________________________________________________________________________
00404 void LCDVToplBase::CompressVerticesList() {
00405   GetVertices()->Compress();
00406 }
00407 //__________________________________________________________________________
00408 
00409 //__________________________________________________________________________
00410 Bool_t LCDVToplBase::CheckKs0Vertex(LCDVToplVRT* vrt) {
00411   Bool_t fk0s=kFALSE;
00412 
00413   Int_t    ntrks     = vrt->GetTrackEntries();
00414   Double_t netCharge = vrt->GetVertexCharge();
00415   if (ntrks == 2 && netCharge == 0.) {
00416     Double_t reconMass = vrt->GetMass();
00417     if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) fk0s=kTRUE;
00418   }
00419   
00420   return fk0s;
00421 }
00422 //__________________________________________________________________________
00423 
00424 //__________________________________________________________________________
00425 Bool_t LCDVToplBase::CheckLambda0Vertex(LCDVToplVRT* vrt) {
00426   Bool_t flambda=kFALSE;
00427 
00428   Double_t mpi = 0.13957;
00429   Double_t mp  = 0.93827;
00430 
00431   Int_t    ntrks     = vrt->GetTrackEntries();
00432   Double_t netCharge = vrt->GetVertexCharge();
00433   if (ntrks == 2 && netCharge == 0.) {
00434     TVector3 p3trk1=vrt->GetTrackMomentumVector(0);
00435     TVector3 p3trk2=vrt->GetTrackMomentumVector(1);
00436     TLorentzVector p4trk1a(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mpi*mpi));
00437     TLorentzVector p4trk2a(p3trk2,TMath::Sqrt(p3trk1.Mag2()+mp *mp ));
00438     TLorentzVector p4vrta=p4trk1a+p4trk2a;
00439     if ( TMath::Abs(p4vrta.M() - 1.1157) < m_SeedVrtLambda0Mass ) {
00440       flambda=kTRUE;
00441     } else {
00442       TLorentzVector p4trk1b(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mp *mp ));
00443       TLorentzVector p4trk2b(p3trk2,TMath::Sqrt(p3trk1.Mag2()+mpi*mpi));
00444       TLorentzVector p4vrtb=p4trk1b+p4trk2b;
00445       if ( TMath::Abs(p4vrtb.M() - 1.1157) < m_SeedVrtLambda0Mass ) 
00446         flambda=kTRUE;
00447     }
00448   }
00449   
00450   return flambda;
00451 }
00452 //__________________________________________________________________________
00453 
00454 //__________________________________________________________________________
00455 Int_t LCDVToplBase::GetVrtTrackEntries(){
00456   return m_nvrttrack;
00457 }
00458 //__________________________________________________________________________
00459 
00460 //__________________________________________________________________________
00461 LCDVToplTRK* LCDVToplBase::GetVrtTrackAt(Int_t itrk) {
00462   LCDVToplTRK* trk=0;
00463   if (itrk < m_nvrttrack) {
00464     trk=(LCDVToplTRK*)GetTracks()->At(m_vrttrack_list[itrk]);
00465   }
00466   return trk;
00467 }
00468 //__________________________________________________________________________
00469 
00470 
00471 //__________________________________________________________________________
00472 void LCDVToplBase::AddVrtTrack(LCDVToplTRK* trk) {
00473   if(m_nvrttrack < 100) {
00474     m_vrttrack_list[m_nvrttrack++]=GetTracks()->IndexOf(trk);
00475   } else {
00476     printf("LCDVToplGhost::AddVrtTrack m_nvrttrack >= 100\n");
00477   }
00478 }
00479 //__________________________________________________________________________
00480 
00481 //__________________________________________________________________________
00482 //Get Vertex Charge
00483 Double_t LCDVToplBase::GetVertexCharge() {
00484   Int_t nvtrks=GetVrtTrackEntries();
00485   Double_t q_vertex=0;
00486   for (Int_t itrk=0 ; itrk < nvtrks ; itrk++) {
00487     q_vertex += GetVrtTrackAt(itrk)->GetCharge();
00488   }
00489   return q_vertex;
00490 }
00491 //__________________________________________________________________________
00492 
00493 
00494 //--------------------------------------------------------------------------
00495 // Heavy flavor jet tag relates...
00496 //--------------------------------------------------------------------------
00497 
00498 //__________________________________________________________________________
00499 Int_t LCDVToplBase::GetSeedVertexByClassicMethod() {
00500   Int_t nseedvrt = -1;
00501   Int_t nvrt = GetVertexEntries();  
00502   if (nvrt < 2) return nseedvrt;
00503 
00504   Double_t maxvsig = -1.;
00505   LCDVToplVRT* vrt=GetVertexAt(0);
00506   for(Int_t ivrt =1; ivrt < nvrt ; ivrt++){
00507     vrt = GetVertexAt(ivrt);    
00508     
00509     // Apply cuts to select seed vertex candidates
00510 
00511     // Cut 1
00512     Double_t decayL  = vrt->GetDecayLength(m_ipv);
00513     if (decayL > m_SeedVrtRMax) continue;
00514 
00515     // Cut 2
00516     Int_t ntrks = vrt->GetTrackEntries();
00517     Double_t netCharge = vrt->GetVertexCharge();
00518     if (ntrks == 2 && netCharge == 0.) {      
00519       Double_t reconMass = vrt->GetMass();
00520       if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) continue;
00521     }
00522 
00523     // Cut 3
00524     TVector3 vrtpos(vrt->GetVertexVector() - m_ipv);
00525     if (vrtpos*m_pja < 0.) continue;
00526 
00527     // Cut 4
00528     TVector3 vrtmon(vrt->GetFourVector().Vect());
00529     if (vrtpos*vrtmon < 0.) continue;
00530 
00531     TMatrixD iper(3,3); GetBeamIPErrorMatrix(iper);
00532     Double_t decayLerr = TMath::Sqrt(vrt->GetErrorDecayLength(m_ipv,iper));
00533     Double_t vrtsig = decayL/decayLerr;
00534     if (vrtsig < m_SeedVrtSigMin) continue;
00535 
00536     // Select the closest vertex
00537     Double_t vsig=vrt->GetVsig();
00538     if (nseedvrt < 0 || vsig > maxvsig) {
00539       maxvsig = vsig;  nseedvrt  = ivrt;
00540     }
00541   }  
00542 
00543   return nseedvrt;
00544 }
00545 //__________________________________________________________________________
00546 
00547 //__________________________________________________________________________
00548 Bool_t LCDVToplBase::GetTrackAttachByClassicMethod(LCDVToplTRK* ltrk,
00549                                                    LCDVToplVRT* seedvrt,
00550                                                    TVector3* xdir) {
00551   Double_t dist = (seedvrt->GetVertexVector() - m_ipv).Mag();
00552 
00553   Double_t l    = 0;
00554   Double_t lodi = 0;
00555   Double_t trdi = 0;
00556   Double_t anta = 0;
00557   TVector3 xseedvrt;
00558   if (xdir) {
00559     xseedvrt=*xdir;
00560   } else {
00561     xseedvrt=seedvrt->GetVertexVector(); xseedvrt += m_ipv;
00562   }
00563 
00564   ltrk->GetTrdiLodiAnta(l,m_ipv,xseedvrt,&trdi,&lodi,&anta);
00565   Double_t ldst = lodi/dist;
00566 
00567   // Apply cuts to select track attachment
00568   Bool_t f_ret=kTRUE;
00569   if (trdi  > m_TrdiMax ) f_ret=kFALSE;
00570   if (lodi  < m_LodiMin ) f_ret=kFALSE;
00571   if (ldst  < m_LodrMin ) f_ret=kFALSE;
00572   if (ldst  > m_LodrMax ) f_ret=kFALSE;
00573   if (anta  > m_AntaMax ) f_ret=kFALSE;
00574   if (trdi  > lodi      ) f_ret=kFALSE;
00575 
00576   return f_ret;
00577 }
00578 //__________________________________________________________________________
00579 
00580 //__________________________________________________________________________
00581 void LCDVToplBase::SetPtrVertexTrack() {
00582   LCDVToplTRK* trk=0;
00583   LCDVToplVRT* vrt=0;
00584 
00585   // Set pointers (track->vertex)
00586   Int_t itrk=0;
00587   Int_t ntrk=m_ntrk2vrt;
00588   for (itrk=0 ; itrk < m_ntrk2vrt ; itrk++) {
00589     trk=GetTrackAt(itrk);
00590     if (m_trk2vrt[itrk] >= 0) {
00591       vrt=GetVertexAt(m_trk2vrt[itrk]);
00592       trk->SetVertexEntries(1);
00593       trk->SetVertexAt(vrt,0);
00594     } else {
00595       trk->SetVertexEntries(0);
00596     }
00597   }
00598   
00599   // Set pointers (vertex->track)
00600   Int_t ivrt=0;
00601   Int_t nvrt=GetVertexEntries();
00602   Int_t jtrk=0;
00603   for (ivrt=0; ivrt < nvrt ; ivrt++) {
00604     vrt =GetVertexAt(ivrt);
00605     ntrk=vrt->GetTrackEntries();
00606     for (itrk=0 ; itrk < ntrk ; itrk++) {
00607       vrt->SetTrackAt(GetTrackAt(m_vrt2trk[jtrk++]),itrk);
00608     }
00609   }
00610   if (jtrk != m_nvrt2trk) {
00611     printf("???? jtrk != m_nvrt2trk in LCDVToplBase::SetPtrVertexTrack\n");
00612     printf("???? jtrk = %d m_nvrt2trk=%d\n",jtrk,m_nvrt2trk);
00613   }
00614 }
00615 //__________________________________________________________________________
00616 
00617 //__________________________________________________________________________
00618 /*
00619 //__________________________________________________________________________
00620 void LCDVToplBase::SortVertex() {
00621   Int_t    nvrt=GetVertexEntries();
00622   Int_t    ivrti,ivrtj;
00623   Double_t disti,distj;
00624   LCDVToplVRT* vrti;
00625   LCDVToplVRT* vrtj;
00626   LCDVToplVRT  vrtk;
00627   for (ivrti=0 ; ivrti < nvrt-1 ; ivrti++) {
00628     for (ivrtj=ivrti+1 ; ivrtj < nvrt ; ivrtj++) {
00629       vrti=GetVertexAt(ivrti);
00630       vrtj=GetVertexAt(ivrtj);
00631       disti=vrti->GetSortVal();
00632       distj=vrtj->GetSortVal();
00633       if (disti > distj) {
00634         vrtk = *vrtj;
00635         *vrtj = *vrti;
00636         *vrti = vrtk;
00637       }
00638     }
00639   }
00640 }
00641 //__________________________________________________________________________
00642 
00643 //__________________________________________________________________________
00644 void LCDVToplBase::AddVertex(LCDVToplVRT* vrt) {
00645   Int_t nvrt=GetVertexEntries();
00646   Int_t ivrt=0;
00647   Bool_t f_no_same_vertex=kTRUE;
00648   for (ivrt=0 ; ivrt < nvrt-1 ; ivrt++) {
00649     if (vrt == GetVertexAt(ivrt)) {
00650       f_no_same_vertex=kFALSE;
00651       break;
00652     }
00653   }
00654 
00655   if (f_no_same_vertex) {
00656     SortVertexByVsig();
00657   } else {
00658     vertex_list->Remove(vrt);
00659   }
00660 }
00661 //__________________________________________________________________________
00662 
00663 //__________________________________________________________________________
00664 void LCDVToplBase::AddVertexIntoTrack(LCDVToplTRK* trk, LCDVToplVRT* vrt) {
00665   Int_t nvrt=trk->GetVertexEntries();
00666   Int_t ivrt=0;
00667   Bool_t f_no_same_vertex=kTRUE;
00668   LCDVToplVRT* tvrt=0;
00669   for (ivrt=0; ivrt < nvrt ; ivrt++) {
00670     tvrt=(LCDVToplVRT*)trk->GetVertexAt(ivrt);
00671     if (tvrt == vrt) { 
00672       f_no_same_vertex=kFALSE; 
00673       break;
00674     }
00675   }
00676   if (f_no_same_vertex) {
00677     trk->AddVertex(vrt);
00678     SortVertexByVsigInTrack(trk);
00679   }  
00680 }
00681 //__________________________________________________________________________
00682 
00683 //__________________________________________________________________________
00684 //Get IP probabirity
00685 Double_t LCDVToplBase::GetIPProbabirity(TVector3 xvtx) {
00686   TMatrixD ip2xvtx(3,1);
00687   for (Int_t i=0 ; i < 3 ; i++) { ip2xvtx(i,0)=xvtx[i] - ipv[i]; }
00688   TMatrixD tmp0(m_iper_inv,TMatrixD::kMult,ip2xvtx);
00689   TMatrixD tmp1(ip2xvtx,TMatrixD::kTransposeMult,tmp0);
00690   Double_t prfu = 0.0;
00691   // protect against floating overflow
00692   if (tmp1(0,0) < 100.0) {
00693     prfu = TMath::Exp(-0.5*tmp1(0,0));
00694   }
00695 
00696   return prfu;
00697 }
00698 //__________________________________________________________________________
00699 
00700 //__________________________________________________________________________
00701 //Get Initial vertexposition ip and track 
00702 TVector3 LCDVToplBase::GetInitialVertexIP2T(Int_t itrk) {
00703   //initial vertex point for given tracks
00704   LCDVToplTRK* trk=GetTrackAt(itrk);
00705   TVector3 xpt=trk->GetPositionVector(0.0);
00706   Double_t stk2=trk->GetErrorDistance(0.0,ipv[0],ipv[1],ipv[2]);
00707   TVector3 dx=xpt-ipv;
00708   Double_t adx=dx.Mag();
00709   TMatrixD dddxT(3,1);
00710   for (Int_t i=0 ; i < 3 ; i++) { dddxT(i,0)=dx[i]/adx; }
00711   TMatrixD tmp0(m_iper,TMatrixD::kMult,dddxT);
00712   TMatrixD tmp1(dddxT,TMatrixD::kTransposeMult,tmp0);
00713   Double_t sip2=tmp1(0,0);
00714   Double_t wtot=1.0/sip2 + 1.0/stk2;
00715   return ipv + (1.0/wtot/stk2)*dx;
00716 
00717 }
00718 //__________________________________________________________________________
00719 
00720 //__________________________________________________________________________
00721 void LCDVToplBase::MergeVertex(LCDVToplVRT* vrti,LCDVToplVRT* vrtj,Int_t f_fit) {
00722   if (vrti == vrtj) { return; }
00723 
00724   LCDVToplTRK* trk=0;
00725   Int_t ntrk=vrtj->GetTrackEntries();
00726   Int_t itrk=0;
00727   for (itrk=0; itrk < ntrk ; itrk++) {
00728     trk =vrtj->GetTrackAt(itrk);
00729     vrti->AddTrack(trk,f_fit);
00730     trk->RemoveVertex(vrtj);
00731     AddVertexIntoTrack(trk,vrti);
00732     trk->CompressVerticesList();
00733   }
00734 
00735   if (f_fit > 1) {
00736     vrti->SetVsig(GetVertexSignificance(vrti->GetVertexVector()));
00737   }
00738 
00739   RemoveVertex(vrtj);
00740 
00741   CompressVerticesList();
00742 }
00743 //__________________________________________________________________________
00744 
00745 //__________________________________________________________________________
00746 void LCDVToplBase::MergeVertex(Int_t i, Int_t j, Int_t f_fit) {
00747   if (i == j) { return; }
00748 
00749   LCDVToplVRT* vrti=GetVertexAt(i);
00750   LCDVToplVRT* vrtj=GetVertexAt(j);
00751 
00752   MergeVertex(vrti,vrtj,f_fit);
00753 }
00754 //__________________________________________________________________________
00755 
00756 //__________________________________________________________________________
00757 void LCDVToplBase::CheckV0Vertex() {
00758   Int_t ntrk=GetTrackEntries();
00759 
00760   LCDVToplTRK* trki=0;
00761   LCDVToplTRK* trkj=0;
00762   Int_t        itrk=0;
00763   Int_t        jtrk=0;
00764 
00765   TVector3 xvinit;
00766   LCDVToplVRT vrt;
00767   for (itrk=1 ; itrk < ntrk ; itrk++) { 
00768     trki=GetTrackAt(itrk);
00769     
00770     for (jtrk=0 ; jtrk < itrk ;jtrk++) { 
00771       trkj=GetTrackAt(jtrk);
00772 
00773       //track-track combination
00774       xvinit=GetInitialVertexT2T(itrk,jtrk);
00775 
00776       vrt.Clean();
00777       vrt.SetVertex(xvinit);
00778       vrt.AddTrack(trki);
00779       vrt.AddTrack(trkj);
00780 
00781       //V0 check...
00782       if (vrt.GetChi2() < 16.0) {
00783         vrt.CalcFourVector();
00784         if (CheckKs0Vertex(&vrt)) {
00785           trki->SetV0flag(1);
00786           trkj->SetV0flag(1);
00787         }
00788         if (CheckLambda0Vertex(&vrt)) {
00789           trki->SetV0flag(2);
00790           trkj->SetV0flag(2);
00791         }
00792       }
00793 
00794     }
00795   }
00796 }
00797 //__________________________________________________________________________
00798 
00799 
00800 //__________________________________________________________________________
00801 Double_t LCDVToplBase::GetHFTagNN() {
00802   // Get Seed Vertex
00803   Int_t iseedvrt = GetSeedVertex();
00804   return GetHFTagNN(iseedvrt);
00805 }
00806 //__________________________________________________________________________
00807 
00808 //__________________________________________________________________________
00809 Double_t LCDVToplBase::GetHFTagNN(Int_t iseed) {
00810   if (iseed < 0 ) { return -1.0; } //no seed vertex...
00811   if (m_Masstag == 0.0) {
00812     CalcPTMass(iseed);
00813   }
00814   LCDVToplVRT* seedvrt = GetVertexAt(iseed);
00815   TVector3     ipseed(seedvrt->GetVertexVector()-ipv);
00816   Double_t dist=ipseed.Mag();
00817   Double_t mvrt=GetMassTag();
00818   Double_t pvrt=15*mvrt-GetP4Vertex().Rho();
00819   Double_t dvrt=TMath::Log10(dist);
00820   Double_t tvrt=GetVrtTrackEntries();
00821   Double_t nvrt=GetVertexEntriesExceptV0();
00822   Float_t nninp[5];
00823   nninp[0]=TMath::Max(0.0,TMath::Min(1.0,mvrt/ 5.0));
00824   nninp[1]=TMath::Max(0.0,TMath::Min(1.0,(pvrt+20.0)/100.0));
00825   nninp[2]=TMath::Max(0.0,TMath::Min(1.0,(dvrt+2.0)/2.5));
00826   nninp[3]=TMath::Max(0.0,TMath::Min(1.0,tvrt/10.0));
00827   nninp[4]=TMath::Max(0.0,TMath::Min(1.0,(nvrt-1.0)/3.0));
00828   Float_t nnout;
00829   LCDNNHFSelect(nninp,&nnout,0);
00830   return nnout;
00831 }
00832 //__________________________________________________________________________
00833 
00834 //__________________________________________________________________________
00835 void LCDVToplBase::PrintOutForDebug() {
00836 
00837   if (dblv > 19) {
00838     Int_t nvrt=GetVertexEntries();
00839     LCDVToplVRT* vrt;
00840     Int_t        nvtrk;
00841     for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
00842       vrt=GetVertexAt(ivrt);
00843       nvtrk=vrt->GetTrackEntries();
00844       cout<<"    "
00845           <<"ivrt="<<ivrt<<" "
00846           <<"vrt="<<vrt<<" "
00847           <<"ntrk="<<nvtrk<<" "
00848           <<"decayL="<<vrt->GetDecayLength(ipv)<<" "
00849           <<"chi2="<<vrt->GetChi2()<<" "
00850           <<"Vsig="<<vrt->GetVsig()<<" "
00851           <<endl;
00852       if (dblv > 24) {
00853         cout<<"      "
00854             <<"vpos="
00855             <<vrt->GetVertexX()<<" "
00856             <<vrt->GetVertexY()<<" "
00857             <<vrt->GetVertexZ()<<" "
00858             <<endl;
00859       }
00860       if (dblv > 29) {
00861         for (Int_t itrk=0 ; itrk < nvtrk ; itrk++)
00862           cout<<"      "
00863               <<"itrk="<<itrk<<" "
00864               <<"trk="<<vrt->GetTrackAt(itrk)<<" "
00865               <<endl;
00866       }
00867     }
00868 
00869     if (dblv > 34) {
00870       LCDVToplTRK* trk;
00871       for (Int_t itrk=0 ; itrk < GetTrackEntries() ; itrk++) {
00872         trk=(LCDVToplTRK*)GetTrackAt(itrk);
00873         cout<<"      "
00874             <<"itrk="<<itrk<<" "
00875             <<"trk="<<trk<<" "
00876             <<endl;
00877         for (Int_t itvrt=0 ; itvrt < trk->GetVertexEntries() ; itvrt++) {
00878           vrt=(LCDVToplVRT*)(trk->GetVertexAt(itvrt));
00879           cout<<"    "
00880               <<"    "<<" "
00881               <<"itvrt="<<itvrt<<" "
00882               <<"vrt="<<vrt<<" "
00883               <<endl;
00884         }
00885       }
00886     }
00887 
00888   }
00889 
00890 }
00891 //__________________________________________________________________________
00892 */
00893 

Generated on Tue Jul 18 18:33:59 2006 for LCDROOT by  doxygen 1.4.6