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

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------------
00002 // $Id: LCDVToplGhost.cxx,v 1.3 2001/10/10 18:15:39 toshi Exp $
00003 // ----------------------------------------------------------------------------
00004 #include "LCDVToplGhost.h"
00005 
00006 #include "TVirtualFitter.h"
00007 #include "TArrayD.h"
00008 
00009 #include <stdio.h>
00010 #include <iostream.h>
00011 
00012 TVirtualFitter *hFitter=0;
00013 extern void LCDVToplGhostAxisFit(Int_t &npar, Double_t *gin, Double_t &f, 
00014                                  Double_t *u, Int_t flag);
00015 
00016 extern "C"{
00017     Int_t LCDNNSeedVertexGhostZpole(Float_t *in, Float_t *out, Int_t init);
00018     Int_t LCDNNSeedVertexOneProngGhostZpole(Float_t *in, Float_t *out, 
00019                                             Int_t init);
00020     Int_t LCDNNTrackAttachGhostZpole(Float_t *in, Float_t *out, Int_t init);
00021     Int_t LCDNNHFSelectGhostZpole(Float_t *in, Float_t *out, Int_t init);
00022     Int_t LCDNNHFSelectOneProngGhostZpole(Float_t *in, Float_t *out, 
00023                                           Int_t init);
00024     Int_t LCDNNHFSelectNoVertexGhostZpole(Float_t *in, Float_t *out, 
00025                                           Int_t init);
00026 
00027     Int_t LCDNNSeedVertexGhostZH(Float_t *in, Float_t *out, Int_t init);
00028     Int_t LCDNNSeedVertexOneProngGhostZH(Float_t *in, Float_t *out, 
00029                                          Int_t init);
00030     Int_t LCDNNTrackAttachGhostZH(Float_t *in, Float_t *out, Int_t init);
00031 
00032 }
00033 
00034 ClassImp(LCDVToplGhost)
00035 
00036 //__________________________________________________________________________
00037 //Constructor
00038 LCDVToplGhost::LCDVToplGhost() {
00039   m_dblevel=0;
00040   Init();
00041 }
00042 //__________________________________________________________________________
00043 
00044 //__________________________________________________________________________
00045 //Constructor
00046 LCDVToplGhost::LCDVToplGhost(Int_t dblv) {
00047   m_dblevel=dblv;
00048   Init();
00049   m_dblevel=dblv;
00050 }
00051 //__________________________________________________________________________
00052 
00053 //__________________________________________________________________________
00054 //Deconstructor
00055 LCDVToplGhost::~LCDVToplGhost() {
00056   CleanBase(1);
00057 }
00058 //__________________________________________________________________________
00059 
00060 //__________________________________________________________________________
00061 void LCDVToplGhost::Init() {
00062   if (m_dblevel > 99) cout<<" <<LCDVToplGhost::Init>> InitBase...";
00063   InitBase(); //Initialize LCDVToplBase.
00064   if (m_dblevel > 99) cout<<"done"<< endl;
00065 
00066   if (m_dblevel > 99) cout<<" <<LCDVToplGhost::Init>> InitControlParameter...";
00067   InitControlParameter(); //Initialize LCDVToplGhost
00068   if (m_dblevel > 99) cout<<"done"<< endl;
00069 
00070   if (m_dblevel > 99) cout<<" <<LCDVToplGhost>> ...";
00071   m_gaxi.SetXYZ(0,0,0);
00072   m_gwid=m_gini;
00073   if (m_dblevel > 99) cout<<"done" << endl;
00074 }
00075 //__________________________________________________________________________
00076 
00077 //__________________________________________________________________________
00078 void LCDVToplGhost::InitControlParameter() {
00079   m_gini =25e-4;
00080   m_gmin =25e-4;
00081   m_cmax =1.0;
00082   m_pcut =0.01;
00083   m_nvreq=0;
00084   m_momf =0;
00085 }
00086 //__________________________________________________________________________
00087 
00088 //__________________________________________________________________________
00089 void LCDVToplGhost::Clean(Int_t iflg) {
00090   CleanBase(iflg);
00091 }
00092 //__________________________________________________________________________
00093 
00094 //__________________________________________________________________________
00095 //Find ghost track.
00096 Int_t LCDVToplGhost::FindGhostTrack() {
00097   if (m_dblevel > 9) 
00098     cout << " <<LCDVToplGhost::FindGhostTrack>> Start ..."<<endl;
00099 
00100   //
00101   // Count # of tracks and create an array for track->vertex mapping.
00102   //
00103   Int_t itrk=0;
00104   m_ntrk2vrt=0;
00105   for (itrk=0 ; itrk < m_nvtrack; itrk++) {
00106     GetTrackAt(itrk)->SetVertexEntries(0);
00107     m_ntrk2vrt++;
00108   }
00109   if (m_trk2vrt) {
00110     delete [] m_trk2vrt;
00111   }
00112   m_trk2vrt = new Int_t[m_ntrk2vrt];
00113   if (m_trkl) {
00114     delete [] m_trkl;
00115   }
00116   m_trkl=new Double_t[m_ntrk2vrt];
00117   for (itrk=0 ; itrk < m_ntrk2vrt ; itrk++) {
00118     m_trkl[itrk]=0;
00119   }
00120 
00121   // N of track check.
00122   if (m_nvtrack < 1) {
00123     return 0;
00124   }
00125 
00126   // Calcurate initial chi**2 wrt ip.
00127   // for (Int_t ii=0 ; ii < m_nvtrack ; ii++) {
00128   // }
00129 
00130   // Initialize ghost track.
00131   m_gaxi=m_pja.Unit();
00132   Double_t thetav=m_gaxi.Theta();
00133   Double_t phiv  =m_gaxi.Phi();
00134   Double_t cstheta=0;
00135   Double_t sntheta=0;
00136   Double_t csphi=0;
00137   Double_t snphi=0;
00138 
00139   LCDVToplTRK *tkghost= new((*m_vtrack_list)[m_nvtrack++]) LCDVToplTRK();
00140   // SetLCDVToplTRKForGhost(tkghost);
00141 
00142   Double_t arglist[100];
00143   hFitter = TVirtualFitter::Fitter(this);
00144   hFitter->Clear();
00145   hFitter->SetFCN(LCDVToplGhostAxisFit);
00146   arglist[0] = -1;
00147   // arglist[0] = 1;
00148   hFitter->ExecuteCommand("SET PRINT", arglist,1);
00149   arglist[0] = 0;
00150   hFitter->ExecuteCommand("SET NOW",   arglist,0);
00151   arglist[0] = 1;
00152   hFitter->ExecuteCommand("SET ERR",arglist,1);
00153   Int_t ipar=0;
00154   Char_t* parname[2]={(Char_t*)"phi",(Char_t*)"theta"};
00155   Double_t par[2];
00156   par[0]=phiv;
00157   par[1]=thetav;
00158   Double_t epar[2]={0.04,0.04};
00159   Double_t parl[2];
00160   Double_t parh[2];
00161   parl[0]=parl[1]=-2*TMath::Pi();
00162   parh[0]=parh[1]=+2*TMath::Pi();
00163 
00164   Int_t itry=0;
00165   for (itry=0 ; itry < 2 ; itry++) {
00166     for (ipar=0 ; ipar < 2 ; ipar++) {
00167       hFitter->SetParameter(ipar,parname[ipar],par[ipar],epar[ipar],
00168                             parl[ipar],parh[ipar]);
00169     }
00170     arglist[0] = 500;
00171     arglist[1] = 1.;
00172     hFitter->ExecuteCommand("MIGRAD", arglist ,2);
00173 
00174     Char_t parName[50];
00175     for (ipar=0 ; ipar < 2 ; ipar++) {
00176       hFitter->GetParameter(ipar,parName,par[ipar],epar[ipar],
00177                             parl[ipar],parh[ipar]);
00178     }
00179     Double_t gwidmax=-100.0;
00180     phiv  =par[0];
00181     thetav=par[1];
00182     csphi  =TMath::Cos(phiv);
00183     snphi  =TMath::Sin(phiv);
00184     cstheta=TMath::Cos(thetav);
00185     sntheta=TMath::Sin(thetav);
00186     m_gaxi.SetXYZ(sntheta*csphi,sntheta*snphi,cstheta);
00187   
00188     // Calc. here to make sure gwidmax, dtrk at best thetav,phiv
00189     CalcChi2ForGhostTrackAxis(&gwidmax);
00190     
00191     // Calculate gwid do all tracks within chi**2 < cmax
00192     if (gwidmax > 0) { gwidmax = TMath::Sqrt(gwidmax); }
00193     if (gwidmax > m_gmin) {
00194       m_gwid = gwidmax; 
00195     } else {
00196       m_gwid = m_gmin;
00197     }
00198   }
00199  
00200   // End...
00201   m_vtrack_list->Remove(tkghost);
00202   m_nvtrack--;
00203 
00204   if (m_dblevel > 9)
00205     cout << " <<LCDVToplGhost::FindGhostTrack>> End..." << endl;
00206   return 1;
00207 }
00208 //__________________________________________________________________________
00209 
00210 //__________________________________________________________________________
00211 // Calc. chi**2 for ghost track axis determination.
00212 Double_t LCDVToplGhost::CalcChi2ForGhostTrackAxis(Double_t* gwidmax) {
00213   if (gwidmax) *gwidmax = -100.0;
00214 
00215   Double_t chi2sum=0;
00216   LCDVToplTRK* tkghost=GetTrackAt(m_nvtrack-1);
00217 
00218   SetLCDVToplTRKForGhost(tkghost);
00219 
00220   // Loop over all tracks
00221   Double_t     prfi =0;
00222   Double_t     prfj =0;
00223   Double_t     chisq=0;
00224   Double_t     de   =0;
00225   LCDVToplTRK* trk  =0;
00226   Double_t     gwidn=0;
00227   TVector3     xvrt;
00228   Int_t        behindIP;
00229   for (Int_t itrk=0 ; itrk < m_nvtrack-1 ; itrk++) {
00230     trk     = GetTrackAt(itrk);
00231     xvrt    = GetInitialVertexT2T(trk,tkghost,m_trkl[itrk],0.0);
00232 
00233     // LCDVToplVRT vrt(xvrt);
00234     // vrt.AddTrack(tkghost);
00235     // vrt.AddTrack(trk);
00236     // xvrt=vrt.GetVertexVector();
00237 
00238     behindIP=m_gaxi*(xvrt - m_ipv) < 0;
00239     if (behindIP) { // behind IP
00240       prfi    = trk->GetTrackChiSqContrib(m_ipv,m_ipv);
00241       prfj    = tkghost->GetTrackChiSqContrib(m_ipv,m_ipv);
00242       chisq   = prfi + prfj;
00243       if (gwidmax) {
00244         m_trkl[itrk]=0.0;
00245       }
00246     } else {
00247       prfi    = trk->GetTrackChiSqContrib(xvrt,m_ipv);
00248       prfj    = tkghost->GetTrackChiSqContrib(xvrt,m_ipv);
00249       chisq   = prfi + prfj;
00250       if (gwidmax) {
00251         de = m_gwid*TMath::Sqrt(prfj);
00252         if (de > 1e-5) {
00253           gwidn = chisq*TMath::Power(m_gwid,4)*(chisq/m_cmax - 1.0)/de/de
00254             + m_gwid*m_gwid;
00255           if (gwidn > *gwidmax) {
00256             // if (gwidn > 0.0 && gwidn < 0.5 && gwidn > *gwidmax) {
00257             *gwidmax = gwidn; 
00258           }
00259         }
00260         m_trkl[itrk]=trk->GetLPOCA2(xvrt);
00261       }
00262     }
00263     chi2sum+= chisq;
00264   }
00265 
00266   // Try soft jet core weight
00267   Double_t ajet = m_pja.Angle(m_gaxi);
00268   chi2sum += TMath::Power(TMath::Power(TMath::Abs(ajet-0.02),0.8)/0.3,2);
00269   // chi2sum += TMath::Power(TMath::Abs(ajet-0.02)/0.1,2);
00270 
00271   return chi2sum;
00272 }
00273 //__________________________________________________________________________
00274 
00275 //__________________________________________________________________________
00276 // Set LCDVToplTRK for ghost axis.
00277 void LCDVToplGhost::SetLCDVToplTRKForGhost(LCDVToplTRK* tkghost) {
00278 
00279   Double_t tk[5];
00280   Double_t pT= 10000.0; // Put high momentum.
00281   Double_t d0= m_ipv.Perp();
00282   Double_t gPt=m_gaxi.Perp();
00283   if (m_ipv.Cross(m_gaxi).Z() > 0) d0 = -d0;
00284   tk[0]= d0;
00285   tk[1]= gPt > 0 ? m_gaxi.Phi() : 0.0 ;
00286   tk[2]=-m_bfield.Z()/333.567/pT;
00287   tk[3]= m_ipv.Z();
00288   tk[4]= gPt > 1e-5 ? m_gaxi.Z()/m_gaxi.Perp() : TMath::Sign(1e5,m_gaxi.Z());
00289   Double_t etk[15];
00290   Double_t gwid2=m_gwid*m_gwid;
00291   for (Int_t i=0 ; i < 15 ; i++) {
00292     etk[i]=0;
00293   }
00294   etk[ 0]=gwid2;
00295   etk[ 2]=1e-16;
00296   etk[ 5]=1e-16;
00297   etk[ 9]=gwid2*(1.0 + tk[4]*tk[4]);
00298   etk[14]=1e-16;
00299 
00300   tkghost->SetUp(10000,tk,etk,1.0,m_bfield.Z(),-1,m_ipv);
00301 
00302 }
00303 //__________________________________________________________________________
00304 
00305 //__________________________________________________________________________
00306 //Kinematic algorithm option, with ghost track
00307 Int_t LCDVToplGhost::FindVertices(Int_t f_v0flag){
00308 
00309   if (m_dblevel > 9) 
00310     cout << " <<LCDVToplGhost::FindVertices>> start" << endl;
00311 
00312   m_gaxi.SetXYZ(0,0,0);
00313   m_gwid=m_gini;
00314 
00315   DeleteVertices(0); // Clean up all vertices.
00316   if (m_dblevel > 9) 
00317     cout << " <<LCDVToplGhost::FindVertices>> pass  1" << endl;
00318 
00319   //
00320   //V0 check...
00321   //
00322   if (f_v0flag) {
00323     CheckV0Vertex();
00324   }
00325 
00326   FindGhostTrack();// determine ghost track direction.
00327   if (m_dblevel > 9) 
00328     cout << " <<LCDVToplGhost::FindVertices>> pass  2" << endl;
00329   
00330   // Setup Ghost track.
00331   LCDVToplTRK *tkghost= new((*m_vtrack_list)[m_nvtrack++]) LCDVToplTRK();
00332   SetLCDVToplTRKForGhost(tkghost);
00333   if (m_dblevel > 9) 
00334     cout << " <<LCDVToplGhost::FindVertices>> pass  3" << endl;
00335 
00336   //
00337   // Calculate longitudinal momentum of each track w.r.t. a ghost track.
00338   //
00339   Int_t itrk=0;
00340   Int_t ntrk=GetTrackEntries();
00341   TArrayD pl2gst(ntrk);
00342   TVector3 ptrk;
00343   for (itrk=0 ; itrk < ntrk-1 ; itrk++) { // Loop over all tracks.
00344     ptrk = GetTrackAt(itrk)->GetMomentumVector(0.0);
00345     pl2gst[itrk]=TMath::Max(ptrk*m_gaxi,1e-5);
00346   }
00347   if (m_dblevel > 9) 
00348     cout << " <<LCDVToplGhost::FindVertices>> pass  4" << endl;
00349 
00350   TClonesArray& va=*GetVertices();  // Vertices Array
00351 
00352   //
00353   // Make initial candidate vertices.
00354   //
00355   m_nvertex=0;
00356   LCDVToplTRK* trki=0;
00357   LCDVToplVRT* vrti=0;
00358   LCDVToplVRT* vrtj=0;
00359   vrtj=new (va[m_nvertex++]) LCDVToplVRT(m_ipv); // IP.
00360   TMatrixD iper(3,3); GetBeamIPErrorMatrix(iper);
00361   vrtj->SetErrorMatrixVertex(iper);
00362   vrtj->AddTrack(tkghost);
00363   for (itrk=0 ; itrk < ntrk-1 ; itrk++) { // Loop over all tracks.
00364     trki = GetTrackAt(itrk);
00365     vrti = new(va[m_nvertex++]) LCDVToplVRT(GetInitialVertexT2T(trki,tkghost));
00366     vrti->AddTrack(tkghost);
00367     vrti->AddTrack(trki);
00368     if ((vrti->GetVertexVector() - m_ipv)*m_gaxi < 0) { // behind IP
00369       // GetVertices()->Remove(vrti);
00370       // m_nvertex--;
00371       RemoveVertex(vrti);
00372       vrtj->AddTrack(trki);
00373     }
00374   }
00375   if (m_dblevel > 9) 
00376     cout << " <<LCDVToplGhost::FindVertices>> pass  5" << endl;
00377 
00378   //
00379   // The loop for ghost-track vertexing.
00380   //
00381   Int_t ivrti=0;
00382   Int_t ivrtj=0;
00383   LCDVToplVRT vrt;
00384   Double_t pmax=-1.0;
00385   Double_t cl =0;
00386   Double_t c2 =0;
00387   Double_t pc2=0;
00388   Int_t    nvtk=0;
00389   Int_t    ivrtv[2]={0,0};
00390   TVector3 xvrt;
00391   while(1) {
00392     pmax=-1.0;
00393     ivrtv[0]=-1;
00394     ivrtv[1]=-1;
00395     for (ivrti=1 ; ivrti < m_nvertex ; ivrti++) { // Loop over all vertices
00396       for (ivrtj=0 ; ivrtj < ivrti ; ivrtj++) {
00397         vrti = GetVertexAt(ivrti);
00398         vrtj = GetVertexAt(ivrtj);
00399         vrt  =*vrti;
00400         vrt +=*vrtj;
00401         // cl=vrt.GetConfLevel();
00402         xvrt = vrt.GetVertexVector();
00403         nvtk = vrt.GetTrackEntries();
00404         c2   = 0;
00405         for (itrk=0 ; itrk < nvtk ; itrk++) {
00406           trki=vrt.GetTrackAt(itrk);
00407           if (trki == tkghost && ivrtj == 0) {//IP vertex and ghost-track
00408             TMatrixD eipi(3,3); GetBeamIPErrorMatrixInv(eipi);
00409             TMatrixD dvm(3,1);
00410             dvm(0,0)=xvrt.X() - m_ipv.X();
00411             dvm(1,0)=xvrt.Y() - m_ipv.Y();
00412             dvm(2,0)=xvrt.Z() - m_ipv.Z();
00413             TMatrixD tmp1(eipi,TMatrixD::kMult,dvm);
00414             TMatrixD tmp2(dvm,TMatrixD::kTransposeMult,tmp1);
00415             pc2 = tmp2(0,0); // Add IP error
00416           } else {
00417             pc2 = trki->GetTrackChiSqContrib(xvrt,m_ipv);
00418             if (trki != tkghost) {
00419               pc2 += m_momf*TMath::Log(pl2gst[GetTrackIdx(trki)]);
00420             }
00421           }
00422           c2 += pc2;
00423         }
00424         cl=TMath::Prob(c2,2*nvtk-3);
00425         if (cl > pmax) {
00426           pmax    =cl;
00427           ivrtv[0]=ivrti;
00428           ivrtv[1]=ivrtj;
00429         }
00430       }
00431     }
00432 
00433     if (m_dblevel > 9) 
00434       cout << " <<LCDVToplGhost::FindVertices>> pmax=" << pmax
00435            << " m_pcut=" << m_pcut
00436            << " ivrtv=" << ivrtv[0] << "," << ivrtv[1]
00437            << " nvrt=" << m_nvertex
00438            << endl;
00439         
00440     if (pmax < m_pcut) break; // 
00441 
00442     vrti  = GetVertexAt(ivrtv[0]);
00443     vrtj  = GetVertexAt(ivrtv[1]);
00444     *vrtj+=*vrti;
00445     RemoveVertex(vrti);
00446   }
00447   if (m_dblevel > 9) 
00448     cout << " <<LCDVToplGhost::FindVertices>> pass  6" << endl;
00449 
00450   //
00451   // Calculate vertex momentum
00452   //
00453   Int_t gstidx=0;
00454   for (ivrti=0 ; ivrti < m_nvertex ; ivrti++) {
00455     vrti=GetVertexAt(ivrti);
00456     vrti->CalcFourVector();
00457     if ((gstidx = vrti->IndexOfTrack(tkghost)) > -1) {
00458       *vrti->GetFourVectorPtr() -= vrti->GetTrack4MomentumVector(gstidx);
00459       vrti->RemoveTrack(tkghost,0);
00460     }
00461   }
00462   if (m_dblevel > 9) 
00463     cout << " <<LCDVToplGhost::FindVertices>> pass  7" << endl;
00464 
00465   //
00466   // Remove Ghost track
00467   //
00468   m_vtrack_list->Remove(tkghost);
00469   m_nvtrack--;
00470   if (m_dblevel > 9) 
00471     cout << " <<LCDVToplGhost::FindVertices>> pass  8" << endl;
00472 
00473   //
00474   // Calculate vertex distance from IP, projected on a ghost-track.
00475   //
00476   Double_t dtrk;
00477   for (ivrti=0 ; ivrti < m_nvertex ; ivrti++) {
00478     vrti=GetVertexAt(ivrti);
00479     dtrk=(vrti->GetVertexVector()-m_ipv)*m_gaxi;
00480     TVector3 vpos=vrti->GetVertexVector();
00481     if (ivrti > 0 && dtrk <= 0) { // behind IP. should be IP vertex.
00482       *vrtj += *vrti;
00483       vrtj->CalcFourVector();
00484       RemoveVertex(vrti);
00485     } else {
00486       if (ivrti == 0) {
00487         vrtj=vrti;
00488         dtrk=0;
00489       }
00490       vrti->SetSortVal(dtrk);
00491     }
00492   }
00493   if (m_dblevel > 9) 
00494     cout << " <<LCDVToplGhost::FindVertices>> pass  9" << endl;
00495 
00496   //
00497   // Check # of tracks in a vertex
00498   //
00499   for (ivrti=1 ; ivrti < m_nvertex ; ivrti++) {
00500     vrti = GetVertexAt(ivrti);
00501     nvtk = vrti->GetTrackEntries();
00502     if (nvtk > 0) continue;
00503     RemoveVertex(vrti);
00504   }
00505   if (m_dblevel > 9) 
00506     cout << " <<LCDVToplGhost::FindVertices>> pass 10" << endl;
00507 
00508   //
00509   // Sort vertex array
00510   //
00511   GetVertices()->Sort();
00512   if (m_dblevel > 9) 
00513     cout << " <<LCDVToplGhost::FindVertices>> pass 11" << endl;
00514 
00515   for (ivrti=1 ; ivrti < m_nvertex ; ivrti++) {
00516     vrti = GetVertexAt(ivrti);
00517     nvtk = vrti->GetTrackEntries();
00518   }
00519 
00520   //
00521   // Make vertex->track map
00522   //
00523   m_nvrt2trk=0;
00524   for (ivrti=0 ; ivrti < m_nvertex ; ivrti++) {
00525     m_nvrt2trk+=GetVertexAt(ivrti)->GetTrackEntries();
00526   }
00527   if (m_vrt2trk) {
00528     delete [] m_vrt2trk;
00529   }
00530   m_vrt2trk=new Int_t[m_nvrt2trk];
00531 
00532   Int_t nvrt2trk=0;
00533   for (ivrti=0 ; ivrti < m_nvertex ; ivrti++) {
00534     vrti=GetVertexAt(ivrti);
00535     nvtk=vrti->GetTrackEntries();
00536     for (itrk=0 ; itrk < nvtk ; itrk++) {
00537       trki=vrti->GetTrackAt(itrk);
00538       if (trki==0) {
00539         continue;
00540       }
00541       trki->SetVertexEntries(1);
00542       trki->SetVertexAt(vrti,0);
00543       m_vrt2trk[nvrt2trk++]=GetTrackIdx(trki);
00544     }
00545   }
00546 
00547   Int_t ntrk2vrt=0;
00548   ntrk2vrt=0;
00549   for (itrk=0 ; itrk < m_nvtrack; itrk++) {
00550     trki=GetTrackAt(itrk);
00551     if (trki->GetVertexEntries() > 0) {
00552       vrti=trki->GetVertexAt(0);
00553       m_trk2vrt[ntrk2vrt++]=GetVertexIdx(vrti);
00554     } else {
00555       m_trk2vrt[ntrk2vrt++]=-1;
00556     }
00557   }
00558   if (m_dblevel > 9) 
00559     cout << " <<LCDVToplGhost::FindVertices>> pass 12" << endl;
00560 
00561   //
00562   // That's it folks.
00563   //
00564   return m_nvertex;
00565 }
00566 //__________________________________________________________________________
00567 
00568 //__________________________________________________________________________
00569 void LCDVToplGhost::CalcPTMass(Int_t iseed) {
00570   m_Rmass    = 0.;
00571   m_PtMass   = 0.;
00572   m_PtMassEr = 0.;
00573   m_Masstag  = 0.;
00574   m_MissPt   = 0.;
00575   Int_t nvtrk=GetVrtTrackEntries();
00576   Int_t ivtrk=0;
00577   for(ivtrk=0 ; ivtrk < nvtrk ; ivtrk++) {
00578     m_vrttrack_list[ivtrk]=0;
00579   }
00580   m_nvrttrack=0;
00581   m_p4Vertex.SetXYZT(0,0,0,0);
00582   
00583   // Check seed vertex index
00584   if (iseed == 0) iseed=GetSeedVertex();
00585   if (iseed < 0) return;
00586 
00587   // Seed vertex.
00588   LCDVToplVRT* seedvrt=GetVertexAt(iseed);
00589 
00590   LCDVToplVRT* vrt  = 0;
00591   LCDVToplTRK* trk = 0;
00592   Int_t nvrt = GetVertexEntries();
00593   Int_t ivrt = 0;
00594   TVector3 ptrk;
00595   Double_t nnone =0;
00596   // const Double_t Mpi=0.13957018;
00597   if (m_UseNN) {
00598 
00599     for (ivrt=1 ; ivrt < nvrt ; ivrt++) {
00600       vrt   = GetVertexAt(ivrt);
00601       nvtrk = vrt->GetTrackEntries();      
00602       if (nvtrk > 1 && GetNNoutputSeedVertex(vrt) < m_SeedVrtNN) {
00603         continue;
00604       } else if (nvtrk == 1 && ivrt != iseed) { // in case of one-prong vertex
00605         nnone=GetNNoutputTrackAttach(vrt->GetTrackAt(0),iseed);
00606         if (nnone < m_TrackAttNN) continue;
00607       }
00608       m_p4Vertex += *vrt->GetFourVectorPtr();
00609       for (ivtrk=0 ; ivtrk < nvtrk ; ivtrk++) {
00610         trk = vrt->GetTrackAt(ivtrk);
00611         if (ivrt == iseed) {
00612           trk->SetTrkVal(1);
00613         } else {
00614           trk->SetTrkVal(2);
00615         }
00616         AddVrtTrack(trk);
00617       }
00618     }
00619 
00620   } else {
00621 
00622     TVector3 xseedvrt(*GetGhostAxis()); xseedvrt += m_ipv;
00623     for (ivrt=1 ; ivrt < nvrt ; ivrt++) {
00624       vrt         = GetVertexAt(ivrt);
00625       m_p4Vertex += *vrt->GetFourVectorPtr();
00626       nvtrk=vrt->GetTrackEntries();
00627       
00628       for (ivtrk=0 ; ivtrk < nvtrk ; ivtrk++) {
00629         trk = vrt->GetTrackAt(ivtrk);
00630         if (ivtrk == iseed) {
00631           trk->SetTrkVal(1);
00632           AddVrtTrack(trk);
00633         } else {
00634           if (GetTrackAttachByClassicMethod(trk,seedvrt,&xseedvrt)) {
00635             trk->SetTrkVal(2);
00636             AddVrtTrack(trk);
00637           }
00638         }
00639       }
00640     }
00641 
00642   }
00643 
00644   // Raw Mass
00645   m_Rmass = m_p4Vertex.Mag(); 
00646 
00647   // Momentum vector of reconstructed Vertex
00648   TVector3 p3Vertex(m_p4Vertex.Vect());
00649 
00650   // vector?
00651   TVector3 ipseed=seedvrt->GetVertexVector() - m_ipv;
00652   TVector3 bv;
00653   bv=(ipseed*m_gaxi)/m_gaxi.Mag()*m_gaxi.Unit();
00654   bv=ipseed;
00655   Double_t dist  =bv.Mag();
00656 
00657   // Missing Pt
00658   Double_t pt = p3Vertex.Perp(bv);
00659   m_MissPt = pt;
00660 
00661   // Pt corrected mass
00662   m_PtMass = TMath::Sqrt(m_Rmass*m_Rmass + pt*pt) + TMath::Abs(pt);
00663 
00664   //Find delta phi and delta lambda.
00665   Double_t sedxy = bv.Perp();
00666   Double_t pxy   = p3Vertex.Pt();
00667   Double_t dPhi  = p3Vertex.DeltaPhi(bv);
00668   Double_t dLam  = TMath::ATan((p3Vertex.Z()/pxy - bv.Z()/sedxy)/
00669                        (1.+p3Vertex.Z()*bv.Z()/(pxy*sedxy)));
00670 
00671   //Rotate vertex error matrix 
00672   TMatrixD evrt(3,3); seedvrt->GetErrorMatrixVertex(evrt);
00673   //TVector3 pvn   = p3Vertex;
00674   TVector3 pvn   = bv;
00675   Double_t pvnphi= pvn.Phi();
00676   Double_t csphi = TMath::Cos(pvnphi);
00677   Double_t snphi = TMath::Sin(pvnphi);
00678   Double_t csthe = pvn.CosTheta();
00679   Double_t snthe = TMath::Sqrt(1.0 - csthe*csthe);
00680   TMatrixD dxpdxT(3,3);
00681   dxpdxT(0,0)=+snphi; dxpdxT(0,1)=-csphi*csthe; dxpdxT(0,2)= csphi*snthe;
00682   dxpdxT(1,0)=-csphi; dxpdxT(1,1)=-snphi*csthe; dxpdxT(1,2)= snphi*snthe;
00683   dxpdxT(2,0)=+0.0  ; dxpdxT(2,1)=       snthe; dxpdxT(2,2)=       csthe;
00684   //dxpdxT(0,0)=+snphi; dxpdxT(1,0)=-csphi*csthe; dxpdxT(2,0)= csphi*snthe;
00685   //dxpdxT(0,1)=-csphi; dxpdxT(1,1)=-snphi*csthe; dxpdxT(2,1)= snphi*snthe;
00686   //dxpdxT(0,2)=+0.0  ; dxpdxT(1,2)=       snthe; dxpdxT(2,2)=       csthe;
00687   TMatrixD tmp0(evrt,TMatrixD::kMult,dxpdxT);
00688   TMatrixD verrot(tmp0,TMatrixD::kTransposeMult,dxpdxT);
00689 
00690   //Take errors into accound
00691   //Allow sigma IP error to make conservative Min missing P
00692   //Also allow seed vertex error n sigma Min missing pt
00693 
00694   Double_t dPhiM = TMath::Abs(dPhi) - 
00695     TMath::Sqrt(m_PTIP[0]*m_PTIP[0] + m_NVsigma*m_NVsigma*verrot(0,0)
00696          + m_NVsigma*m_NVsigma*verrot(2,2)
00697                 *TMath::Sin(dPhi)*TMath::Sin(dPhi) )/sedxy;
00698   if (dPhiM < 0.) dPhiM = 0.;
00699 
00700   Double_t dLamM = TMath::Abs(dLam) - 
00701     TMath::Sqrt(m_PTIP[1]*m_PTIP[1]*sedxy*sedxy/(dist*dist) 
00702                 + m_NVsigma*m_NVsigma*verrot(1,1)
00703                 + m_NVsigma*m_NVsigma
00704                 *verrot(2,2)*TMath::Sin(dLam)*TMath::Sin(dLam))/dist;
00705   if (dLamM < 0.) dLamM = 0.;
00706 
00707   pt = TMath::Sqrt(TMath::Power(p3Vertex.Mag()*TMath::Sin(dLamM),2)
00708                    +TMath::Power(pxy*TMath::Sin(dPhiM),2));
00709   
00710   // Ptcorrected mass with IP errors
00711   m_PtMassEr = TMath::Sqrt(m_Rmass*m_Rmass + pt*pt) + TMath::Abs(pt); 
00712 
00713   // Pt corrected IP err with 2*raw mass protection
00714   if (m_PtMassEr > 2.*m_Rmass) {
00715     m_Masstag = 2.*m_Rmass; 
00716   } else {
00717     m_Masstag = m_PtMassEr; 
00718   }
00719 
00720 }
00721 //__________________________________________________________________________
00722 
00723 //__________________________________________________________________________
00724 Double_t LCDVToplGhost::GetNNoutputSeedVertex(Int_t ivrt) {
00725   return GetNNoutputSeedVertex(GetVertexAt(ivrt));
00726 }
00727 //__________________________________________________________________________
00728 
00729 //__________________________________________________________________________
00730 Double_t LCDVToplGhost::GetNNoutputSeedVertex(LCDVToplVRT* vrt) {
00731     
00732   // Apply cuts to select seed vertex candidates
00733   Double_t decayL  = vrt->GetDecayLength(m_ipv);
00734   TVector3 vrtpos(vrt->GetVertexVector() - m_ipv);
00735   TVector3 vrtmon(vrt->GetFourVector().Vect());
00736   TMatrixD iper(3,3); GetBeamIPErrorMatrix(iper);
00737   Double_t decayLerr = TMath::Sqrt(vrt->GetErrorDecayLength(m_ipv,iper));
00738   Double_t vrtsig = decayL/decayLerr;
00739   Double_t pdang  = vrtpos.Angle(vrtmon);
00740   Int_t    ntrks  = vrt->GetTrackEntries();
00741   Double_t rmvr   = vrt->GetMass();  
00742   decayL=TMath::Log10(decayL);
00743   vrtsig=TMath::Log10(vrtsig);
00744   pdang =TMath::Log10(TMath::Max(pdang,1e-8));
00745   Float_t nninp[5];
00746   Float_t nnout;
00747   if (ntrks > 1) {
00748     nninp[0]=TMath::Min(TMath::Max((decayL+3.0)/ 3.5,0.0),1.0);
00749     nninp[1]=TMath::Min(TMath::Max((vrtsig+0.5)/ 3.5,0.0),1.0);
00750     nninp[2]=TMath::Min(TMath::Max((pdang +3.5)/ 4.0,0.0),1.0);
00751     nninp[3]=TMath::Min(TMath::Max((rmvr  +0.0)/ 7.0,0.0),1.0);
00752     nninp[4]=TMath::Min(TMath::Max((ntrks +0.0)/10.0,0.0),1.0);
00753     // LCDNNSeedVertexGhost(nninp,&nnout,0);
00754     switch (m_UseNN) {
00755     case 1: // Z-pole
00756       LCDNNSeedVertexGhostZpole(nninp,&nnout,0);
00757       break;
00758     case 2: // ZH @ 500GeV
00759       LCDNNSeedVertexGhostZH(nninp,&nnout,0);
00760       break;
00761     default:
00762       break;
00763     }
00764   } else {
00765     nninp[0]=TMath::Min(TMath::Max((decayL+3.0)/ 3.5,0.0),1.0);
00766     nninp[1]=TMath::Min(TMath::Max((vrtsig+0.5)/ 3.5,0.0),1.0);
00767     nninp[2]=TMath::Min(TMath::Max((pdang +3.5)/ 4.0,0.0),1.0);
00768     switch (m_UseNN) {
00769     case 1: // Z-pole
00770       LCDNNSeedVertexOneProngGhostZpole(nninp,&nnout,0);
00771       break;
00772     case 2: // ZH @ 500GeV
00773       LCDNNSeedVertexOneProngGhostZH(nninp,&nnout,0);
00774       break;
00775     default:
00776       break;
00777     }
00778   }
00779 
00780   return nnout;
00781 }
00782 //__________________________________________________________________________
00783 
00784 //__________________________________________________________________________
00785 Int_t LCDVToplGhost::GetSeedVertex() {
00786   Int_t seedvertex_id=0;
00787   switch(m_UseNN) {
00788   case 0:
00789     seedvertex_id=GetSeedVertexByClassicMethod();
00790     break;
00791   default:
00792     seedvertex_id=GetSeedVertexByNNMethod();
00793     break;
00794   }
00795   
00796   return seedvertex_id;
00797 }
00798 //__________________________________________________________________________
00799 
00800 //__________________________________________________________________________
00801 Int_t LCDVToplGhost::GetSeedVertexByNNMethod() {
00802   Int_t nseedvrt = -1;
00803   Int_t nvrt = GetVertexEntries();  
00804   if (nvrt < 2) return nseedvrt;
00805 
00806   LCDVToplVRT* vrt=0;
00807   Double_t maxnn = -1.;
00808   Int_t    ntrks =0;
00809   Float_t  nnout;
00810   Int_t    ivrt=0;
00811   for(ivrt =1; ivrt < nvrt ; ivrt++){
00812 
00813     vrt = GetVertexAt(ivrt);
00814     ntrks=vrt->GetTrackEntries();
00815     if (ntrks == 1) continue; // skip one-prong vertex
00816 
00817     // Cut 1
00818     if (CheckKs0Vertex(vrt)) continue;
00819     // if (CheckLambda0Vertex(vrt)) continue;
00820 
00821     nnout=GetNNoutputSeedVertex(ivrt);
00822     if (nnout < m_SeedVrtNN) continue;
00823 
00824     // Select
00825     if (nseedvrt < 0 || nnout > maxnn) {
00826       maxnn = nnout;  nseedvrt  = ivrt;
00827     }
00828 
00829   }
00830 
00831   if (nseedvrt < 0) { // Then look one-prong vertex
00832     for(ivrt=0 ; ivrt < nvrt ; ivrt++) {
00833       vrt   = GetVertexAt(ivrt);
00834       ntrks = vrt->GetTrackEntries();
00835       if (ntrks != 1) continue; // skip non-one-prong vertex
00836       if (vrt->GetTrackAt(0)->GetV0flag()) continue; // Check V0 track
00837       
00838       nnout=GetNNoutputSeedVertex(ivrt);
00839       if (nnout < m_SeedVrt1PrngNN) continue;
00840 
00841       // Select
00842       if (nseedvrt < 0 || nnout > maxnn) {
00843         maxnn = nnout;  nseedvrt  = ivrt;
00844       }
00845 
00846     }
00847   }
00848 
00849   return nseedvrt;
00850 
00851 }
00852 //__________________________________________________________________________
00853 
00854 //__________________________________________________________________________
00855 Double_t LCDVToplGhost::GetNNoutputTrackAttach(LCDVToplTRK* ltrk, 
00856                                                Int_t iseedvrt) {
00857   Double_t dist = (GetVertexAt(iseedvrt)->GetVertexVector() - m_ipv).Mag();
00858   
00859   Double_t l    = 0;
00860   Double_t lodi = 0;
00861   Double_t trdi = 0;
00862   Double_t anta = 0;
00863   TVector3 xseedvrt(*GetGhostAxis()); xseedvrt+=m_ipv;
00864   ltrk->GetTrdiLodiAnta(l,m_ipv,xseedvrt,&trdi,&lodi,&anta);
00865   Double_t ldst = lodi/dist;
00866   TMatrixD iper(3,3); GetBeamIPErrorMatrix(iper);
00867   Double_t sbds = 
00868     TMath::Log10(ltrk->GetDistanceNormByError(l,m_ipv,iper));
00869 
00870   //Apply cuts to select track attachment
00871   Float_t nninp[6];
00872   nninp[0]=TMath::Max(0.0,TMath::Min(trdi/0.3        ,1.0));
00873   nninp[1]=TMath::Max(0.0,TMath::Min((lodi-0.3)/3.0  ,1.0));
00874   nninp[2]=TMath::Max(0.0,TMath::Min((ldst-0.5)/3.0  ,1.0));
00875   nninp[3]=TMath::Max(0.0,TMath::Min(anta/TMath::Pi(),1.0));
00876   nninp[4]=TMath::Max(0.0,TMath::Min(sbds/3.0        ,1.0));
00877   Float_t nnout;
00878   switch (m_UseNN) {
00879   case 1: //Z-pole
00880     LCDNNTrackAttachGhostZpole(nninp,&nnout,0);
00881     break;
00882   case 2: //ZH
00883     LCDNNTrackAttachGhostZH(nninp,&nnout,0);
00884     break;
00885   }
00886 
00887   return nnout;
00888 }
00889 //__________________________________________________________________________
00890 
00891 //__________________________________________________________________________
00892 void LCDVToplGhost::CalcHFTagNN(Int_t iseed) {
00893   m_HFBTagNN=-1.0;
00894   m_HFCTagNN=-1.0;
00895   m_HF1Prong=-1;
00896   if (iseed < 0 ) { //no seed vertex...
00897     CalcHFTagNN_NoVertex(); 
00898     return;
00899   }
00900   if (iseed == 0 || m_Masstag == 0.0) {
00901     CalcPTMass(iseed);
00902   }
00903   if(iseed <= 0) {
00904     CalcHFTagNN_NoVertex();
00905     return;
00906   }
00907 
00908   Int_t ivrt=0;
00909 
00910   LCDVToplVRT* seedvrt = GetVertexAt(iseed);
00911   TVector3     ipseed(seedvrt->GetVertexVector()-m_ipv);
00912   Double_t dist=ipseed.Mag();
00913   Double_t mvrt=GetMassTag();
00914   Double_t pvrt=15*mvrt-GetP4Vertex().Rho();
00915   Double_t dvrt=TMath::Log10(dist);
00916   Int_t    tvrt=GetVrtTrackEntries();
00917   TMatrixD eipv(3,3); GetBeamIPErrorMatrix(eipv);
00918   Double_t edln = TMath::Sqrt(seedvrt->GetErrorDecayLength(m_ipv,eipv));
00919   Double_t bist = dist/edln;
00920   Double_t bvrt = TMath::Log10(bist);
00921   
00922   //  Double_t nvrt=GetVertexEntriesExceptV0();
00923   Int_t nvrt=1;
00924   Int_t novt=0;
00925   Int_t ntrk=0;
00926   Int_t nvrt0=GetVertexEntries();
00927   Double_t nnvrt=0;
00928   LCDVToplVRT* vrt=0;
00929   for(ivrt=1 ; ivrt < nvrt0 ; ivrt++) {
00930     if (ivrt == iseed) continue;
00931     vrt=GetVertexAt(ivrt);
00932     ntrk=vrt->GetTrackEntries();
00933     if (ntrk > 1) {
00934       nnvrt=GetNNoutputSeedVertex(ivrt);
00935       if (nnvrt < m_SeedVrtNN) continue;
00936     } else {
00937       nnvrt=GetNNoutputTrackAttach(vrt->GetTrackAt(0),iseed);
00938       if (nnvrt < m_TrackAttNN) continue;
00939     }
00940     nvrt++;
00941     if (ntrk == 1) novt++;
00942   }
00943   Float_t nninp[7];
00944   Float_t nnout[2];
00945   if (tvrt > 1) {
00946     nninp[0]=TMath::Max(0.0,TMath::Min(1.0,mvrt/ 7.0));
00947     nninp[1]=TMath::Max(0.0,TMath::Min(1.0,(pvrt+20.0)/80.0));
00948     nninp[2]=TMath::Max(0.0,TMath::Min(1.0,(dvrt+2.0)/2.5));
00949     nninp[3]=TMath::Max(0.0,TMath::Min(1.0,tvrt/10.0));
00950     nninp[4]=TMath::Max(0.0,TMath::Min(1.0,(nvrt-1.0)/4.0));
00951     nninp[5]=TMath::Max(0.0,TMath::Min(1.0,(novt-0.0)/5.0));
00952     nninp[6]=TMath::Max(0.0,TMath::Min(1.0,(bvrt-0.5)/2.5));
00953     switch (m_UseNN) {
00954     case 1: // Z-pole
00955       LCDNNHFSelectGhostZpole(nninp,nnout,0);
00956       break;
00957     case 2: // ZH @ 500 GeV
00958       break;
00959     default:
00960       break;
00961     }
00962     m_HFBTagNN=nnout[0];
00963     m_HFCTagNN=nnout[1];
00964     m_HF1Prong=0;
00965   } else if (tvrt == 1) {
00966     LCDVToplTRK* ltrk=seedvrt->GetTrackAt(0);
00967     Double_t btrk = TMath::Log10(ltrk->GetDistanceNormByError(0.0,m_ipv,eipv));
00968     Double_t ptrk = ltrk->GetMomentumVector(0.0).Mag()/m_ptk.Mag();
00969     TVector3 xseedvrt(*GetGhostAxis());
00970     xseedvrt+= m_ipv;
00971     Double_t l    = 0;
00972     Double_t lodi = 0;
00973     Double_t trdi = 0;
00974     Double_t anta = 0;
00975     ltrk->GetTrdiLodiAnta(l,m_ipv,xseedvrt,&trdi,&lodi,&anta);
00976     anta=TMath::Log10(TMath::Max(anta,1e-2));
00977     nninp[0]=TMath::Max(0.0,TMath::Min(1.0,(dvrt+2.0)/2.5));
00978     nninp[1]=TMath::Max(0.0,TMath::Min(1.0,(bvrt-0.5)/3.0));
00979     nninp[2]=TMath::Max(0.0,TMath::Min(1.0,(btrk-0.5)/3.0));
00980     nninp[3]=TMath::Max(0.0,TMath::Min(1.0,(anta+2.0)/2.5));
00981     nninp[4]=TMath::Max(0.0,TMath::Min(1.0,(ptrk-0.0)/0.5));
00982     switch (m_UseNN) {
00983     case 1: // Z-pole
00984       LCDNNHFSelectOneProngGhostZpole(nninp,nnout,0);
00985       break;
00986     case 2: // ZH @ 500 GeV
00987       break;
00988     default:
00989       break;
00990     }
00991     m_HFBTagNN=nnout[0];
00992     m_HFCTagNN=nnout[1];
00993     m_HF1Prong=1;
00994   } else {
00995     printf(" ....Error something is wrong in LCDVToplGhost::CalcHFTagNN\n");
00996     printf(" tvrt = %d\n",tvrt);
00997     m_HFBTagNN=-1;
00998     m_HFCTagNN=-1;
00999     m_HF1Prong=0;
01000   }
01001   return;
01002 }
01003 //__________________________________________________________________________
01004 
01005 
01006 //__________________________________________________________________________
01007 void LCDVToplGhost::Doit() {
01008   if (m_dblevel > 0) {
01009     printf(" <<LCDVToplGhost::Doit>> FindVertices...start\n");
01010   }
01011   FindVertices(1);             // find vertices with v0 track check.
01012   if (m_dblevel > 0) {
01013     printf(" <<LCDVToplGhost::Doit>> FindVertices...done\n");
01014   }
01015 
01016   if (m_dblevel > 0) {
01017     printf(" <<LCDVToplGhost::Doit>> GetSeedVertex...start\n");
01018   }
01019   Int_t iseed=GetSeedVertex(); // find seed vertex.
01020   if (m_dblevel > 0) {
01021     printf(" <<LCDVToplGhost::Doit>> GetSeedVertex...done iseed=%d nvrt=%d\n",
01022            iseed,GetVertexEntries());
01023   }
01024 
01025   if (m_dblevel > 0) {
01026     printf(" <<LCDVToplGhost::Doit>> CalcPTMass...start\n");
01027   }
01028   CalcPTMass(iseed);           // Calculate PT corrected vertex mass
01029   if (m_dblevel > 0) {
01030     printf(" <<LCDVToplGhost::Doit>> CalcPTMass...done masstag=%f\n",
01031            GetMassTag());
01032   }
01033 
01034   if (m_dblevel > 0) {
01035     printf(" <<LCDVToplGhost::Doit>> CalcHFTagNN...start\n");
01036   }
01037   CalcHFTagNN(iseed);          // Calculate NN outputs for HF tag
01038   if (m_dblevel > 0) {
01039     printf(" <<LCDVToplGhost::Doit>> CalcHFTagNN...done nnb=%f nnc=%f\n",
01040            GetHFBTagNN(),GetHFCTagNN());
01041   }
01042 
01043 }
01044 //__________________________________________________________________________
01045 
01046 //__________________________________________________________________________
01047 Int_t LCDVToplGhost::GetVertexEntriesWithNN(Int_t iseed){
01048   Int_t nvrt=1;
01049   Int_t ntrk=0;
01050   Int_t nvrt0=GetVertexEntries();
01051   LCDVToplVRT* vrt=0;
01052   for(Int_t ivrt=1 ; ivrt < nvrt0 ; ivrt++) {
01053     if (ivrt == iseed) continue;
01054     vrt=GetVertexAt(ivrt);
01055     ntrk=vrt->GetTrackEntries();
01056     if (ntrk > 1) {
01057       if (GetNNoutputSeedVertex(ivrt) < m_SeedVrtNN) continue;
01058     } else {
01059       if (iseed) {
01060         if (GetNNoutputTrackAttach(vrt->GetTrackAt(0),iseed) < m_TrackAttNN) {
01061           continue;
01062         }
01063       } else {
01064         if (GetNNoutputSeedVertex(ivrt) < m_SeedVrt1PrngNN) continue;
01065       }
01066     }
01067     nvrt++;
01068   }
01069   return nvrt;
01070 }
01071 //__________________________________________________________________________
01072 
01073 //__________________________________________________________________________
01074 // Calc. chi**2 for ghost track axis determination.
01075 void LCDVToplGhostAxisFit(Int_t &npar, Double_t *gin, Double_t &f, 
01076                           Double_t *u, Int_t flag) {
01077 
01078   LCDVToplGhost* a=(LCDVToplGhost*)hFitter->GetObjectFit();
01079 
01080   Double_t phiv   =u[0];
01081   Double_t thetav =u[1];
01082   Double_t csphi  =TMath::Cos(phiv);
01083   Double_t snphi  =TMath::Sin(phiv);
01084   Double_t cstheta=TMath::Cos(thetav);
01085   Double_t sntheta=TMath::Sin(thetav);
01086   a->GetGhostAxis()->SetXYZ(sntheta*csphi,sntheta*snphi,cstheta);
01087   f=a->CalcChi2ForGhostTrackAxis(0);
01088 }
01089 //__________________________________________________________________________
01090 
01091 
01092 //==========================================================================
01093 //
01094 //
01095 //
01096 // Below is copied from LCDVToplBase.
01097 //
01098 //
01099 //
01100 //==========================================================================
01101 
01102 
01103 //__________________________________________________________________________
01104 void LCDVToplGhost::SetUp(const TVector3& ipv, 
01105                          const TVector3& bfldv) {
01106 
01107   if (m_dblevel > 9) 
01108     cout << " <<LCDVToplGhost::SetUp1>> Start..." << endl;
01109 
01110   //This loop to make sure position and momentum in trk at 2d POCA
01111   // and  to make LCDVToplTRK.
01112   TVector3 pja(0.0,0.0,0.0); 
01113   //  if (m_pja.Mag2() == 0 && m_nvtrack ) {
01114   if (m_nvtrack) {
01115     for (Int_t i=0; i < m_nvtrack ; i++) {
01116       //Accumurate track momentum.
01117       pja   += GetTrackAt(i)->GetMomentumVector(0.0);
01118     } //loop over zxtrks
01119   }
01120 
01121   SetUp(ipv,bfldv,pja);
01122 
01123   if (m_dblevel > 9) 
01124     cout << " <<LCDVToplGhost::SetUp1>> End..." << endl;
01125 }
01126 //__________________________________________________________________________
01127 
01128 //__________________________________________________________________________
01129 void LCDVToplGhost::SetUp(const TVector3& ipv, 
01130                          const TVector3& bfldv,
01131                          const TVector3& jeta) {
01132   if (m_dblevel > 9) 
01133     cout << " <<LCDVToplGhost::SetUp2>> Start..." << endl;
01134 
01135   CleanBase(0);
01136 
01137   SetIP(ipv);
01138   SetBField(bfldv);
01139   m_pja=jeta;
01140   CalcIPErrorMatrix();
01141 
01142   //This loop to make sure position and momentum in trk at 2d POCA
01143   // and  to make LCDVToplTRK.
01144   m_ptk.SetXYZ(0,0,0);
01145   if (m_nvtrack) {
01146     for (Int_t i=0; i < m_nvtrack ; i++) {
01147       //Accumurate track momentum.
01148       m_ptk += GetTrackAt(i)->GetMomentumVector(0.0);
01149     } //loop over zxtrks
01150   }
01151 
01152   if (m_dblevel > 9) 
01153     cout << " <<LCDVToplGhost::SetUp2>> End..." << endl;
01154 }
01155 //__________________________________________________________________________
01156 
01157 //__________________________________________________________________________
01158 //Init...
01159 void LCDVToplGhost::InitBase() {
01160   m_nvtrack       = 0;
01161   m_vertex_list   = new TClonesArray("LCDVToplVRT",100);
01162   m_nvertex       = 0;
01163   m_vtrack_list   = new TClonesArray("LCDVToplTRK",100);
01164   m_nvrttrack     = 0;
01165   for(Int_t i=0 ; i < 100 ; i++) {
01166     m_vrttrack_list[i]=0;
01167   }
01168     
01169   m_ripe  = 3.0e-4;
01170   m_zipe  = 6.0e-4;
01171 
01172   // Initialize jet total momentum
01173   m_pja.SetXYZ(0.0,0.0,0.0); 
01174   m_ptk.SetXYZ(0.0,0.0,0.0);
01175 
01176   //Set Beam Position (Temporary)  
01177   m_ipv.SetXYZ(0.0,0.0,0.0); 
01178 
01179   // m_dblevel=0;
01180 
01181   //
01182   // Pt corrected mass relates.
01183   //
01184   m_Rmass    = 0.;
01185   m_PtMass   = 0.;
01186   m_PtMassEr = 0.;
01187   m_Masstag  = 0.;
01188   m_MissPt   = 0.;
01189 
01190   m_SeedVrtRMax  =2.3  ;
01191   m_SeedVrtK0Mass=0.005;
01192   m_SeedVrtLambda0Mass=0.005;
01193   m_SeedVrtSigMin=2.0  ;
01194 
01195   m_TrdiMax=0.1;
01196   m_LodiMin=0.05;
01197   m_LodrMin=0.25;
01198   m_LodrMax=100.0;
01199   m_AntaMax=TMath::Pi()/2.0;
01200 
01201   m_PTIP[0]=1.5*m_ripe;
01202   m_PTIP[1]=1.5*m_zipe;
01203   m_NVsigma=1.0;
01204 
01205   m_UseNN          =  0;
01206   m_SeedVrtNN      =  0.3;
01207   m_SeedVrt1PrngNN =  0.7;
01208   m_TrackAttNN     =  0.5;
01209   m_SelectBNN      =  0.5;
01210   m_SelectCNN      =  0.5;
01211   m_HFBTagNN       =-10.0; // NN output for b-tag.
01212   m_HFCTagNN       =-10.0; // NN output for c-tag.
01213   m_HF1Prong       =- 1  ; // =1: one-prong tag/ =0: other
01214 
01215   m_p4Vertex.SetXYZT(0.,0.,0.,0.);
01216 
01217   m_mcjetflav=0;
01218 
01219   m_ntrk2vrt=0;
01220   m_trk2vrt=0;
01221   m_nvrt2trk=0;
01222   m_vrt2trk=0;
01223   m_trkl=0;
01224 
01225 }
01226 //__________________________________________________________________________
01227 
01228 //__________________________________________________________________________
01229 void LCDVToplGhost::DeleteTracks(Int_t iflg){
01230   if (m_vtrack_list)   m_vtrack_list->Delete();
01231   m_nvtrack=0;
01232   if (iflg) {
01233     delete m_vtrack_list;
01234     m_vtrack_list=0;
01235   }
01236   m_nvrttrack=0;
01237   for(Int_t i=0 ; i < 100 ; i++) {
01238     m_vrttrack_list[i]=0;
01239   }
01240 }
01241 //__________________________________________________________________________
01242 
01243 //__________________________________________________________________________
01244 void LCDVToplGhost::DeleteVertices(Int_t iflg){
01245   if (m_vertex_list) m_vertex_list->Delete();
01246   if (iflg) {
01247     delete m_vertex_list;
01248     m_vertex_list=0;
01249   }
01250 }
01251 //__________________________________________________________________________
01252 
01253 //__________________________________________________________________________
01254 //Delete LCDVToplVRT, LCDVToplTRK...
01255 void LCDVToplGhost::CleanBase(Int_t iflg) {
01256   DeleteTracks(iflg);
01257   DeleteVertices(iflg);
01258   if (m_trk2vrt) {
01259     delete [] m_trk2vrt;
01260     m_trk2vrt=0;
01261   }
01262   if (m_vrt2trk) {
01263     delete [] m_vrt2trk;
01264     m_vrt2trk=0;
01265   }
01266   if (m_trkl) {
01267     delete [] m_trkl;
01268     m_trkl=0;
01269   }
01270   m_ntrk2vrt=0;
01271   m_nvrt2trk=0;
01272 }
01273 //__________________________________________________________________________
01274 
01275 //__________________________________________________________________________
01276 void LCDVToplGhost::CalcIPErrorMatrix(){
01277   Double_t ripe2=m_ripe*m_ripe;
01278   Double_t zipe2=m_zipe*m_zipe;
01279   m_iper[0]=ripe2/2.0;
01280   m_iper[1]=0.0      ; m_iper[2]=ripe2/2.0;
01281   m_iper[3]=0.0      ; m_iper[4]=0.0      ; m_iper[5]=zipe2;
01282   TMatrixD t(3,3);
01283   GetBeamIPErrorMatrix(t);
01284   Double_t det_iper;
01285   t.Invert(&det_iper);
01286   SetBeamIPErrorMatrixInv(t);
01287 }
01288 //__________________________________________________________________________
01289 
01290 //__________________________________________________________________________
01291 void LCDVToplGhost::GetBeamIPErrorMatrix(TMatrixD& t) {
01292   t(0,0)=m_iper[0]; t(0,1)=m_iper[1]; t(0,2)=m_iper[3];
01293   t(1,0)=m_iper[1]; t(1,1)=m_iper[2]; t(1,2)=m_iper[4];
01294   t(2,0)=m_iper[3]; t(2,1)=m_iper[4]; t(2,2)=m_iper[5];
01295 }
01296 //__________________________________________________________________________
01297 
01298 //__________________________________________________________________________
01299 void LCDVToplGhost::GetBeamIPErrorMatrixInv(TMatrixD& t) {
01300   t(0,0)=m_iper_inv[0]; t(0,1)=m_iper_inv[1]; t(0,2)=m_iper_inv[3];
01301   t(1,0)=m_iper_inv[1]; t(1,1)=m_iper_inv[2]; t(1,2)=m_iper_inv[4];
01302   t(2,0)=m_iper_inv[3]; t(2,1)=m_iper_inv[4]; t(2,2)=m_iper_inv[5];
01303 }
01304 //__________________________________________________________________________
01305 
01306 //__________________________________________________________________________
01307 void LCDVToplGhost::SetBeamIPErrorMatrix(TMatrixD& t) {
01308   m_iper[0]=t(0,0);
01309   m_iper[1]=t(1,0); m_iper[2]=t(1,1);
01310   m_iper[3]=t(2,0); m_iper[4]=t(2,1); m_iper[5]=t(2,2);
01311 }
01312 //__________________________________________________________________________
01313 
01314 //__________________________________________________________________________
01315 void LCDVToplGhost::SetBeamIPErrorMatrixInv(TMatrixD& t) {
01316   m_iper_inv[0]=t(0,0);
01317   m_iper_inv[1]=t(1,0); m_iper_inv[2]=t(1,1);
01318   m_iper_inv[3]=t(2,0); m_iper_inv[4]=t(2,1); m_iper_inv[5]=t(2,2);
01319 }
01320 //__________________________________________________________________________
01321 
01322 //__________________________________________________________________________
01323 LCDVToplTRK* LCDVToplGhost::GetTrackAt(Int_t itrk) {
01324   return (LCDVToplTRK*)(GetTracks()->At(itrk));
01325 }
01326 //__________________________________________________________________________
01327 
01328 //__________________________________________________________________________
01329 Int_t LCDVToplGhost::GetTrackIdx(LCDVToplTRK* trk) {
01330   return GetTracks()->IndexOf(trk);
01331 }
01332 //__________________________________________________________________________
01333 
01334 //__________________________________________________________________________
01335 //Get # of tracks
01336 Int_t LCDVToplGhost::GetTrackEntries(){
01337   return m_nvtrack;
01338   // return GetTracks()->GetEntries();
01339 }
01340 //__________________________________________________________________________
01341 
01342 //__________________________________________________________________________
01343 LCDVToplVRT* LCDVToplGhost::GetVertexAt(Int_t ivrt) {
01344   return (LCDVToplVRT*)(GetVertices()->At(ivrt));
01345 }
01346 //__________________________________________________________________________
01347 
01348 //__________________________________________________________________________
01349 Int_t LCDVToplGhost::GetVertexIdx(LCDVToplVRT* vrt) {
01350   return GetVertices()->IndexOf(vrt);
01351 }
01352 //__________________________________________________________________________
01353 
01354 //__________________________________________________________________________
01355 Int_t LCDVToplGhost::GetVertexEntries(){
01356   return m_nvertex;
01357   //    return GetVertices()->GetEntries();
01358 }
01359 //__________________________________________________________________________
01360 
01361 //__________________________________________________________________________
01362 Int_t LCDVToplGhost::GetVertexEntriesExceptV0(){
01363   Int_t nvrt=GetVertexEntries();
01364   Int_t nvrtnov0=0;
01365   LCDVToplVRT* vrt;
01366   for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
01367     vrt=GetVertexAt(ivrt);
01368     if (CheckKs0Vertex(vrt))     continue;
01369     if (CheckLambda0Vertex(vrt)) continue;
01370     nvrtnov0++;
01371   }
01372   return nvrtnov0;
01373 }
01374 //__________________________________________________________________________
01375 
01376 //__________________________________________________________________________
01377 LCDVToplTRK* LCDVToplGhost::AddTrack(Int_t tkid, 
01378                                     Double_t* tkpar, 
01379                                     Double_t* e_tkpar, 
01380                                     Double_t chrg, 
01381                                     Double_t mfld, 
01382                                     Int_t mc_ind, 
01383                                     const TVector3& ipv) {
01384   LCDVToplTRK* trk = new((*m_vtrack_list)[m_nvtrack++]) 
01385     LCDVToplTRK(tkid,tkpar,e_tkpar,chrg,mfld,mc_ind,ipv);
01386   //  m_pja += trk->GetMomentumVector(0.0);
01387   m_ptk += trk->GetMomentumVector(0.0);
01388   return trk;
01389 }
01390 //__________________________________________________________________________
01391 
01392 //__________________________________________________________________________
01393 //Get Initial vertexposition rrack and track 
01394 TVector3 LCDVToplGhost::GetInitialVertexT2T(LCDVToplTRK* trk1, 
01395                                            LCDVToplTRK* trk2,
01396                                            Double_t l1, Double_t l2) {
01397   //initial vertex point for given tracks
01398 
01399   TVector3 posv1   = trk1->GetPositionVector(l1);
01400   TVector3 momv1   = trk1->GetMomentumVector(l1);
01401   TVector3 ip2pos1(m_ipv - posv1);
01402   if (l1 != 0.0) {
01403     Double_t mu1   = momv1.X()*ip2pos1.X() + momv1.Y()*ip2pos1.Y();
01404     mu1 /= (momv1.X()*momv1.X() + momv1.Y()*momv1.Y());
01405     posv1 += mu1*momv1;
01406   }
01407   // Double_t omega1  = trk1->GetTrackParameter(2);
01408   // Double_t phi0pr1 = trk1->GetTrackParameter(1) - TMath::Pi()/2.0 + l1/omega1;
01409   Double_t phi0pr1 = momv1.Phi() - TMath::Pi()/2.0;
01410   Double_t cspr1   = TMath::Cos(phi0pr1);
01411   Double_t snpr1   = TMath::Sin(phi0pr1);
01412   Double_t tanlpr1 = trk1->GetTrackParameter(4);
01413   Double_t z0pr1   = posv1.Z();
01414   Double_t x0pr1   = cspr1*(posv1[0] - m_ipv[0]) + snpr1*(posv1[1] - m_ipv[1]);
01415   Double_t s1spr1  = trk1->GetErrorTsiApprox(l1);
01416   Double_t s2spr1  = trk1->GetErrorEtaApprox(l1)*(1.0 + tanlpr1*tanlpr1);
01417 
01418   TVector3 posv2   = trk2->GetPositionVector(l2);
01419   TVector3 momv2   = trk2->GetMomentumVector(l2);
01420   TVector3 ip2pos2(m_ipv - posv2);
01421   if (l2 != 0.0) {
01422     Double_t mu2   = momv2.X()*ip2pos2.X() + momv2.Y()*ip2pos2.Y();
01423     mu2 /= (momv2.X()*momv2.X() + momv2.Y()*momv2.Y());
01424     posv2 += mu2*momv2;
01425   }
01426   // Double_t omega2  = trk2->GetTrackParameter(2);
01427   // Double_t phi0pr2 = trk2->GetTrackParameter(1) - TMath::Pi()/2.0 + l2/omega2;
01428   Double_t phi0pr2 = momv2.Phi() - TMath::Pi()/2.0;
01429   Double_t cspr2   = TMath::Cos(phi0pr2);
01430   Double_t snpr2   = TMath::Sin(phi0pr2);
01431   Double_t tanlpr2 = trk2->GetTrackParameter(4);
01432   Double_t z0pr2   = posv2.Z();
01433   Double_t x0pr2   = cspr2*(posv2[0] - m_ipv[0]) + snpr2*(posv2[1] - m_ipv[1]);
01434   Double_t s1spr2  = trk2->GetErrorTsiApprox(l2);
01435   Double_t s2spr2  = trk2->GetErrorEtaApprox(l2)*(1.0 + tanlpr2*tanlpr2);
01436 
01437   Double_t c1 = 
01438     cspr1*cspr1/s1spr1 + tanlpr1*tanlpr1*snpr1*snpr1/s2spr1 +
01439     cspr2*cspr2/s1spr2 + tanlpr2*tanlpr2*snpr2*snpr2/s2spr2;
01440   Double_t c2 =
01441     snpr1*snpr1/s1spr1 + tanlpr1*tanlpr1*cspr1*cspr1/s2spr1 +
01442     snpr2*snpr2/s1spr2 + tanlpr2*tanlpr2*cspr2*cspr2/s2spr2;
01443   Double_t c3 = 1/s2spr1 + 1/s2spr2;
01444   Double_t c4 =
01445     cspr1*snpr1/s1spr1 - tanlpr1*tanlpr1*snpr1*cspr1/s2spr1 +
01446     cspr2*snpr2/s1spr2 - tanlpr2*tanlpr2*snpr2*cspr2/s2spr2;
01447   Double_t c5 = tanlpr1*snpr1/s2spr1 + tanlpr2*snpr2/s2spr2;
01448   Double_t c6 =-(tanlpr1*cspr1/s2spr1 + tanlpr2*cspr2/s2spr2);
01449   Double_t c7 =
01450     cspr1*x0pr1/s1spr1 + z0pr1*tanlpr1*snpr1/s2spr1 +
01451     cspr2*x0pr2/s1spr2 + z0pr2*tanlpr2*snpr2/s2spr2;
01452   Double_t c8 =
01453     snpr1*x0pr1/s1spr1 - z0pr1*tanlpr1*cspr1/s2spr1 +
01454     snpr2*x0pr2/s1spr2 - z0pr2*tanlpr2*cspr2/s2spr2;
01455   Double_t c9 =z0pr1/s2spr1 + z0pr2/s2spr2;
01456 
01457   TMatrixD rot(3,3);
01458   rot(0,0) = c1; rot(0,1) = c4; rot(0,2)= c5;
01459   rot(1,0) = c4; rot(1,1) = c2; rot(1,2)= c6;
01460   rot(2,0) = c5; rot(2,1) = c6; rot(2,2)= c3;
01461   TMatrixD tmp(3,1);
01462   tmp(0,0) = c7;
01463   tmp(1,0) = c8;
01464   tmp(2,0) = c9;
01465   
01466   TVector3 xdv=m_ipv;
01467   Double_t det_rot;
01468   rot.Invert(&det_rot);
01469   if (det_rot) {
01470     TMatrixD tmp1(rot,TMatrixD::kMult,tmp);
01471     // tmp1.Mult(rot,tmp);
01472     xdv.SetXYZ(tmp1(0,0)+m_ipv.X(),tmp1(1,0)+m_ipv.Y(),tmp1(2,0));
01473     // xdv.SetXYZ(tmp1(0,0)+m_ipv[0],tmp1(1,0)+m_ipv[1],tmp1(2,0));
01474     // avoid being too far off
01475     // TVector3 xdv2(xdv[0],xdv[1],xdv[2]-m_ipv[2]);
01476     // if ( xdv2.Mag2() > 16.0 ) {
01477     //  xdv=m_ipv;
01478     // }
01479   }
01480 
01481   return xdv;
01482 }
01483 //__________________________________________________________________________
01484 
01485 //__________________________________________________________________________
01486 void LCDVToplGhost::RemoveVertexAt(Int_t i) {
01487   RemoveVertex(GetVertexAt(i));
01488 }
01489 //__________________________________________________________________________
01490 
01491 //__________________________________________________________________________
01492 void LCDVToplGhost::RemoveVertex(LCDVToplVRT* vrt) {
01493   GetVertices()->Remove(vrt);
01494   CompressVerticesList();
01495   m_nvertex--;
01496 }
01497 //__________________________________________________________________________
01498 
01499 //__________________________________________________________________________
01500 void LCDVToplGhost::CompressVerticesList() {
01501   GetVertices()->Compress();
01502 }
01503 //__________________________________________________________________________
01504 
01505 //__________________________________________________________________________
01506 Bool_t LCDVToplGhost::CheckKs0Vertex(LCDVToplVRT* vrt) {
01507   Bool_t fk0s=kFALSE;
01508 
01509   Int_t    ntrks     = vrt->GetTrackEntries();
01510   Double_t netCharge = vrt->GetVertexCharge();
01511   if (ntrks == 2 && netCharge == 0.) {
01512     Double_t reconMass = vrt->GetMass();
01513     if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) fk0s=kTRUE;
01514   }
01515   
01516   return fk0s;
01517 }
01518 //__________________________________________________________________________
01519 
01520 //__________________________________________________________________________
01521 Bool_t LCDVToplGhost::CheckLambda0Vertex(LCDVToplVRT* vrt) {
01522   Bool_t flambda=kFALSE;
01523 
01524   Double_t mpi = 0.13957;
01525   Double_t mp  = 0.93827;
01526 
01527   Int_t    ntrks     = vrt->GetTrackEntries();
01528   Double_t netCharge = vrt->GetVertexCharge();
01529   if (ntrks == 2 && netCharge == 0.) {
01530     TVector3 p3trk1=vrt->GetTrackMomentumVector(0);
01531     TVector3 p3trk2=vrt->GetTrackMomentumVector(1);
01532     TLorentzVector p4trk1a(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mpi*mpi));
01533     TLorentzVector p4trk2a(p3trk2,TMath::Sqrt(p3trk2.Mag2()+mp *mp ));
01534     TLorentzVector p4vrta=p4trk1a+p4trk2a;
01535     if ( TMath::Abs(p4vrta.M() - 1.1157) < m_SeedVrtLambda0Mass ) {
01536       flambda=kTRUE;
01537     } else {
01538       TLorentzVector p4trk1b(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mp *mp ));
01539       TLorentzVector p4trk2b(p3trk2,TMath::Sqrt(p3trk2.Mag2()+mpi*mpi));
01540       TLorentzVector p4vrtb=p4trk1b+p4trk2b;
01541       if ( TMath::Abs(p4vrtb.M() - 1.1157) < m_SeedVrtLambda0Mass ) 
01542         flambda=kTRUE;
01543     }
01544   }
01545   
01546   return flambda;
01547 }
01548 //__________________________________________________________________________
01549 
01550 //__________________________________________________________________________
01551 Int_t LCDVToplGhost::GetVrtTrackEntries(){
01552     return m_nvrttrack;
01553 }
01554 //__________________________________________________________________________
01555 
01556 //__________________________________________________________________________
01557 LCDVToplTRK* LCDVToplGhost::GetVrtTrackAt(Int_t itrk) {
01558   LCDVToplTRK* trk=0;
01559   if (itrk < m_nvrttrack) {
01560     trk=(LCDVToplTRK*)GetTracks()->At(m_vrttrack_list[itrk]);
01561   }
01562   return trk;
01563 }
01564 //__________________________________________________________________________
01565 
01566 
01567 //__________________________________________________________________________
01568 void LCDVToplGhost::AddVrtTrack(LCDVToplTRK* trk) {
01569   if(m_nvrttrack < 100) {
01570     m_vrttrack_list[m_nvrttrack++]=GetTracks()->IndexOf(trk);
01571   } else {
01572     printf("LCDVToplGhost::AddVrtTrack m_nvrttrack >= 100\n");
01573   }
01574 }
01575 //__________________________________________________________________________
01576 
01577 //__________________________________________________________________________
01578 //Get Vertex Charge
01579 Double_t LCDVToplGhost::GetVertexCharge() {
01580   Int_t nvtrks=GetVrtTrackEntries();
01581   Double_t q_vertex=0;
01582   for (Int_t itrk=0 ; itrk < nvtrks ; itrk++) {
01583     q_vertex += GetVrtTrackAt(itrk)->GetCharge();
01584   }
01585   return q_vertex;
01586 }
01587 //__________________________________________________________________________
01588 
01589 
01590 //--------------------------------------------------------------------------
01591 // Heavy flavor jet tag relates...
01592 //--------------------------------------------------------------------------
01593 
01594 //__________________________________________________________________________
01595 Int_t LCDVToplGhost::GetSeedVertexByClassicMethod() {
01596   Int_t nseedvrt = -1;
01597   Int_t nvrt = GetVertexEntries();  
01598   if (nvrt < 2) return nseedvrt;
01599 
01600   Double_t maxvsig = -1.;
01601   LCDVToplVRT* vrt=GetVertexAt(0);
01602   for(Int_t ivrt =1; ivrt < nvrt ; ivrt++){
01603     vrt = GetVertexAt(ivrt);    
01604     
01605     // Apply cuts to select seed vertex candidates
01606 
01607     // Cut 1
01608     Double_t decayL  = vrt->GetDecayLength(m_ipv);
01609     if (decayL > m_SeedVrtRMax) continue;
01610 
01611     // Cut 2
01612     Int_t ntrks = vrt->GetTrackEntries();
01613     Double_t netCharge = vrt->GetVertexCharge();
01614     if (ntrks == 2 && netCharge == 0.) {      
01615       Double_t reconMass = vrt->GetMass();
01616       if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) continue;
01617     }
01618 
01619     // Cut 3
01620     TVector3 vrtpos(vrt->GetVertexVector() - m_ipv);
01621     if (vrtpos*m_pja < 0.) continue;
01622 
01623     // Cut 4
01624     TVector3 vrtmon(vrt->GetFourVector().Vect());
01625     if (vrtpos*vrtmon < 0.) continue;
01626 
01627     TMatrixD iper(3,3); GetBeamIPErrorMatrix(iper);
01628     Double_t decayLerr = TMath::Sqrt(vrt->GetErrorDecayLength(m_ipv,iper));
01629     Double_t vrtsig = decayL/decayLerr;
01630     if (vrtsig < m_SeedVrtSigMin) continue;
01631 
01632     // Select the closest vertex
01633     Double_t vsig=vrt->GetVsig();
01634     if (nseedvrt < 0 || vsig > maxvsig) {
01635       maxvsig = vsig;  nseedvrt  = ivrt;
01636     }
01637   }  
01638 
01639   return nseedvrt;
01640 }
01641 //__________________________________________________________________________
01642 
01643 //__________________________________________________________________________
01644 Bool_t LCDVToplGhost::GetTrackAttachByClassicMethod(LCDVToplTRK* ltrk,
01645                                                    LCDVToplVRT* seedvrt,
01646                                                    TVector3* xdir) {
01647   Double_t dist = (seedvrt->GetVertexVector() - m_ipv).Mag();
01648 
01649   Double_t l    = 0;
01650   Double_t lodi = 0;
01651   Double_t trdi = 0;
01652   Double_t anta = 0;
01653   TVector3 xseedvrt;
01654   if (xdir) {
01655     xseedvrt=*xdir;
01656   } else {
01657     xseedvrt=seedvrt->GetVertexVector(); xseedvrt += m_ipv;
01658   }
01659 
01660   ltrk->GetTrdiLodiAnta(l,m_ipv,xseedvrt,&trdi,&lodi,&anta);
01661   Double_t ldst = lodi/dist;
01662 
01663   // Apply cuts to select track attachment
01664   Bool_t f_ret=kTRUE;
01665   if (trdi  > m_TrdiMax ) f_ret=kFALSE;
01666   if (lodi  < m_LodiMin ) f_ret=kFALSE;
01667   if (ldst  < m_LodrMin ) f_ret=kFALSE;
01668   if (ldst  > m_LodrMax ) f_ret=kFALSE;
01669   if (anta  > m_AntaMax ) f_ret=kFALSE;
01670   if (trdi  > lodi      ) f_ret=kFALSE;
01671 
01672   return f_ret;
01673 }
01674 //__________________________________________________________________________
01675 
01676 //__________________________________________________________________________
01677 void LCDVToplGhost::SetPtrVertexTrack() {
01678   LCDVToplTRK* trk=0;
01679   LCDVToplVRT* vrt=0;
01680 
01681   // Set pointers (track->vertex)
01682   Int_t itrk=0;
01683   Int_t ntrk=m_ntrk2vrt;
01684   for (itrk=0 ; itrk < m_ntrk2vrt ; itrk++) {
01685     trk=GetTrackAt(itrk);
01686     if (m_trk2vrt[itrk] >= 0) {
01687       vrt=GetVertexAt(m_trk2vrt[itrk]);
01688       trk->SetVertexEntries(1);
01689       trk->SetVertexAt(vrt,0);
01690     } else {
01691       trk->SetVertexEntries(0);
01692     }
01693   }
01694   
01695   // Set pointers (vertex->track)
01696   Int_t ivrt=0;
01697   Int_t nvrt=GetVertexEntries();
01698   Int_t jtrk=0;
01699   for (ivrt=0; ivrt < nvrt ; ivrt++) {
01700     vrt =GetVertexAt(ivrt);
01701     ntrk=vrt->GetTrackEntries();
01702     for (itrk=0 ; itrk < ntrk ; itrk++) {
01703       vrt->SetTrackAt(GetTrackAt(m_vrt2trk[jtrk++]),itrk);
01704     }
01705   }
01706   if (jtrk != m_nvrt2trk) {
01707     printf("???? jtrk != m_nvrt2trk in LCDVToplGhost::SetPtrVertexTrack\n");
01708     printf("???? jtrk = %d m_nvrt2trk=%d\n",jtrk,m_nvrt2trk);
01709   }
01710 }
01711 //__________________________________________________________________________
01712 
01713 //__________________________________________________________________________
01714 void LCDVToplGhost::CheckV0Vertex() {
01715   Int_t ntrk=GetTrackEntries();
01716 
01717   LCDVToplTRK* trki=0;
01718   LCDVToplTRK* trkj=0;
01719   Int_t        itrk=0;
01720   Int_t        jtrk=0;
01721 
01722   TVector3 xvinit;
01723   LCDVToplVRT vrt;
01724   for (itrk=1 ; itrk < ntrk ; itrk++) { 
01725     trki=GetTrackAt(itrk);    
01726     for (jtrk=0 ; jtrk < itrk ;jtrk++) { 
01727       trkj=GetTrackAt(jtrk);
01728       if (trki->GetCharge()*trkj->GetCharge() >= 0) continue;
01729 
01730       //track-track combination
01731       xvinit=GetInitialVertexT2T(trki,trkj);
01732 
01733       vrt.Clean();
01734       vrt.SetVertex(xvinit);
01735       vrt.AddTrack(trki);
01736       vrt.AddTrack(trkj);
01737 
01738       //V0 check...
01739       if (vrt.GetChi2() < 16.0) {
01740         vrt.CalcFourVector();
01741         if (CheckKs0Vertex(&vrt)) {
01742           trki->SetV0flag(1);
01743           trkj->SetV0flag(1);
01744         }
01745         if (CheckLambda0Vertex(&vrt)) {
01746           trki->SetV0flag(2);
01747           trkj->SetV0flag(2);
01748         }
01749       }
01750 
01751     }
01752   }
01753 }
01754 //__________________________________________________________________________
01755 
01756 //__________________________________________________________________________
01757 void LCDVToplGhost::CalcHFTagNN_NoVertex() {
01758   // To calculate PJ (r-phi), d0/sigma(d0) (1,2), and z0/sigma(z0) (1,2)
01759   LCDVToplTRK* ltrk=0;
01760   TVector3 xtrk;
01761   TVector3 ptrk;
01762   Double_t drtrk,dztrk;
01763   Double_t ertrk,eztrk;
01764   Double_t brtrk,bztrk;
01765   Double_t artrk,aztrk;
01766   Double_t yr=1,yz=1;
01767   Double_t br1mx=0,br2mx=0;
01768   Double_t bz1mx=0,bz2mx=0;
01769   Double_t ar1mx=0,ar2mx=0;
01770   Double_t az1mx=0,az2mx=0;
01771   Double_t ptottk=0;
01772   Int_t itrk=0;
01773   Int_t ntrk1=0,ntrk2=0;
01774   TMatrixD eipv(3,3); GetBeamIPErrorMatrix(eipv);
01775   Double_t l=0;
01776   for (itrk=0 ; itrk < GetTrackEntries() ; itrk++) {
01777     ltrk = GetTrackAt(itrk);
01778     if (ltrk->GetV0flag()) continue; // Skip V0 track
01779     l=ltrk->GetLPOCA2(m_ipv);
01780     xtrk = ltrk->GetPositionVector(l) - m_ipv;
01781     ptrk = ltrk->GetMomentumVector(l);
01782     ptottk=ptrk.Mag();
01783     drtrk = ltrk->GetSignedDistance2D(l,m_ipv,*GetPja());
01784     dztrk = TMath::Abs(xtrk.Z());
01785     if ((*GetPja())*xtrk < 0) {
01786       dztrk = -dztrk;
01787     }
01788     ertrk = ltrk->GetErrorDistance2D(l,m_ipv.X(),m_ipv.Y(),eipv);
01789     eztrk = ltrk->GetErrorDistanceZ(l,m_ipv.Z(),eipv);
01790     brtrk = 0;
01791     if (ertrk > 0) {
01792       brtrk = drtrk/TMath::Sqrt(ertrk);          
01793       if (brtrk > 0) {
01794         if (ntrk1 == 0) {
01795           yr  = TMath::Prob(brtrk*brtrk,1);
01796         } else {
01797           yr *= TMath::Prob(brtrk*brtrk,1);
01798         }
01799         ntrk1++;
01800       }
01801     }
01802     bztrk = 0;
01803     if (eztrk > 0) {
01804       bztrk = dztrk/TMath::Sqrt(eztrk);
01805       if (bztrk > 0) {
01806         if (ntrk2 == 0) {
01807           yz  = TMath::Prob(bztrk*bztrk,1);
01808         } else {
01809           yz *= TMath::Prob(bztrk*bztrk,1);
01810         }
01811         ntrk2++;
01812       }
01813     }
01814     if (ptottk > 2.0) {
01815       if (ertrk > 0) {
01816         artrk=TMath::Abs(brtrk);
01817         if (artrk > ar1mx) {
01818           br2mx = br1mx;
01819           br1mx = brtrk;
01820           ar2mx = ar1mx;
01821           ar1mx = TMath::Abs(brtrk);
01822         } else if (artrk > ar2mx) {
01823           br2mx = brtrk;
01824           ar2mx = TMath::Abs(artrk);
01825         }
01826       }
01827       if (eztrk > 0) {
01828         aztrk=TMath::Abs(bztrk);
01829         if (aztrk > az1mx) {
01830           bz2mx = bz1mx;
01831           bz1mx = bztrk;
01832           az2mx = az1mx;
01833           az1mx = TMath::Abs(bztrk);
01834         } else if (aztrk > az2mx) {
01835           bz2mx = bztrk;
01836           az2mx = TMath::Abs(aztrk);
01837         }
01838       }
01839     }
01840   }
01841   Double_t lnyr=-TMath::Log(TMath::Max(1e-5,yr));
01842   Double_t pp1 =1;
01843   Double_t pj_rf=0;
01844   for (itrk=0 ; itrk < ntrk1 ; itrk++) {
01845     pj_rf += pp1;
01846     pp1   *= lnyr/(itrk+1);
01847   }
01848   pj_rf *= yr;
01849   Double_t lnyz=-TMath::Log(TMath::Max(1e-5,yz));
01850   Double_t pp2 =1;
01851   Double_t pj_rz=0;
01852   for (itrk=0 ; itrk < ntrk2 ; itrk++) {
01853     pj_rz += pp2;
01854     pp2   *= lnyz/(itrk+1);
01855   }
01856   pj_rz *= yz;
01857   Float_t nninp[6];
01858   Float_t nnout[2];
01859   nninp[0]=TMath::Max(0.0,TMath::Min(1.0,pj_rf));
01860   nninp[1]=TMath::Max(0.0,TMath::Min(1.0,pj_rz));
01861   nninp[2]=TMath::Max(0.0,TMath::Min(1.0,(br1mx+10.0)/30.0));
01862   nninp[3]=TMath::Max(0.0,TMath::Min(1.0,(br2mx+10.0)/30.0));
01863   nninp[4]=TMath::Max(0.0,TMath::Min(1.0,(bz1mx+10.0)/30.0));
01864   nninp[5]=TMath::Max(0.0,TMath::Min(1.0,(bz2mx+10.0)/30.0));
01865   
01866   nnout[0]=0;
01867   nnout[1]=0;
01868 
01869   switch (m_UseNN) {
01870   case 1: // Z-pole
01871     LCDNNHFSelectNoVertexGhostZpole(nninp,nnout,0);
01872     break;
01873   case 2: // ZH @ 500 GeV
01874     break;
01875   default:
01876     break;
01877   }
01878 
01879   m_HFBTagNN=nnout[0];
01880   m_HFCTagNN=nnout[1];
01881   m_HF1Prong=2;
01882 }
01883 //__________________________________________________________________________
01884 
01885 //__________________________________________________________________________
01886 /*
01887 //__________________________________________________________________________
01888 //Get IP probabirity
01889 Double_t LCDVToplGhost::GetIPProbabirity(TVector3 xvtx) {
01890   TMatrixD ip2xvtx(3,1);
01891   for (Int_t i=0 ; i < 3 ; i++) { ip2xvtx(i,0)=xvtx[i] - ipv[i]; }
01892   TMatrixD tmp0(m_iper_inv,TMatrixD::kMult,ip2xvtx);
01893   TMatrixD tmp1(ip2xvtx,TMatrixD::kTransposeMult,tmp0);
01894   Double_t prfu = 0.0;
01895   // protect against floating overflow
01896   if (tmp1(0,0) < 100.0) {
01897     prfu = TMath::Exp(-0.5*tmp1(0,0));
01898   }
01899 
01900   return prfu;
01901 }
01902 //__________________________________________________________________________
01903 
01904 //__________________________________________________________________________
01905 //Get Initial vertexposition ip and track 
01906 TVector3 LCDVToplGhost::GetInitialVertexIP2T(Int_t itrk) {
01907   //initial vertex point for given tracks
01908   LCDVToplTRK* trk=GetTrackAt(itrk);
01909   TVector3 xpt=trk->GetPositionVector(0.0);
01910   Double_t stk2=trk->GetErrorDistance(0.0,ipv[0],ipv[1],ipv[2]);
01911   TVector3 dx=xpt-ipv;
01912   Double_t adx=dx.Mag();
01913   TMatrixD dddxT(3,1);
01914   for (Int_t i=0 ; i < 3 ; i++) { dddxT(i,0)=dx[i]/adx; }
01915   TMatrixD tmp0(m_iper,TMatrixD::kMult,dddxT);
01916   TMatrixD tmp1(dddxT,TMatrixD::kTransposeMult,tmp0);
01917   Double_t sip2=tmp1(0,0);
01918   Double_t wtot=1.0/sip2 + 1.0/stk2;
01919   return ipv + (1.0/wtot/stk2)*dx;
01920 
01921 }
01922 //__________________________________________________________________________
01923 
01924 //__________________________________________________________________________
01925 void LCDVToplGhost::PrintOutForDebug() {
01926 
01927   if (dblv > 19) {
01928     Int_t nvrt=GetVertexEntries();
01929     LCDVToplVRT* vrt;
01930     Int_t        nvtrk;
01931     for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
01932       vrt=GetVertexAt(ivrt);
01933       nvtrk=vrt->GetTrackEntries();
01934       cout<<"    "
01935           <<"ivrt="<<ivrt<<" "
01936           <<"vrt="<<vrt<<" "
01937           <<"ntrk="<<nvtrk<<" "
01938           <<"decayL="<<vrt->GetDecayLength(ipv)<<" "
01939           <<"chi2="<<vrt->GetChi2()<<" "
01940           <<"Vsig="<<vrt->GetVsig()<<" "
01941           <<endl;
01942       if (dblv > 24) {
01943         cout<<"      "
01944             <<"vpos="
01945             <<vrt->GetVertexX()<<" "
01946             <<vrt->GetVertexY()<<" "
01947             <<vrt->GetVertexZ()<<" "
01948             <<endl;
01949       }
01950       if (dblv > 29) {
01951         for (Int_t itrk=0 ; itrk < nvtrk ; itrk++)
01952           cout<<"      "
01953               <<"itrk="<<itrk<<" "
01954               <<"trk="<<vrt->GetTrackAt(itrk)<<" "
01955               <<endl;
01956       }
01957     }
01958 
01959     if (dblv > 34) {
01960       LCDVToplTRK* trk;
01961       for (Int_t itrk=0 ; itrk < GetTrackEntries() ; itrk++) {
01962         trk=(LCDVToplTRK*)GetTrackAt(itrk);
01963         cout<<"      "
01964             <<"itrk="<<itrk<<" "
01965             <<"trk="<<trk<<" "
01966             <<endl;
01967         for (Int_t itvrt=0 ; itvrt < trk->GetVertexEntries() ; itvrt++) {
01968           vrt=(LCDVToplVRT*)(trk->GetVertexAt(itvrt));
01969           cout<<"    "
01970               <<"    "<<" "
01971               <<"itvrt="<<itvrt<<" "
01972               <<"vrt="<<vrt<<" "
01973               <<endl;
01974         }
01975       }
01976     }
01977 
01978   }
01979 
01980 }
01981 //__________________________________________________________________________
01982 */
01983 

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