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

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------------
00002 // $Id: LCDVTopl.cxx,v 1.5 2001/10/08 17:53:52 toshi Exp $
00003 // ----------------------------------------------------------------------------
00004 //
00005 ////////////////////////////////////////////////////////////
00006 //
00007 // LCDVTopl
00008 //
00009 //
00010 #include "LCDVTopl.h"
00011 
00012 #include <stdio.h>
00013 #include <iostream.h>
00014 
00015 ClassImp(LCDVTopl)
00016 
00017 //__________________________________________________________________________
00018 //Constructor and Deconstructor
00019 LCDVTopl::LCDVTopl() {
00020   Init();
00021   CalcIPErrorMatrix();
00022 }
00023 //__________________________________________________________________________
00024 
00025 //__________________________________________________________________________
00026 LCDVTopl::LCDVTopl(TObjArray* track_list,
00027                    Int_t ntrk, 
00028                    Int_t* track_indecies,
00029                    const TVector3& ip,
00030                    const TVector3& bfld) {
00031 
00032   Init();
00033   SetUpLCDVTopl(track_list,ntrk,track_indecies,ip,bfld);
00034 }
00035 //__________________________________________________________________________
00036 
00037 //__________________________________________________________________________
00038 LCDVTopl::LCDVTopl(TObjArray* track_list,
00039                    Int_t ntrk, 
00040                    Int_t* track_indecies,
00041                    const TVector3& ip,
00042                    const TVector3& bfld,
00043                    const TVector3& pjet) {
00044   Init();
00045   SetUpLCDVTopl(track_list,ntrk,track_indecies,ip,bfld,pjet);
00046 }
00047 //__________________________________________________________________________
00048 
00049 //__________________________________________________________________________
00050 LCDVTopl::LCDVTopl(TObjArray* track_list,
00051                    const TVector3& ip,
00052                    const TVector3& bfld,
00053                    const TVector3& pjet) {
00054   Init();
00055 
00056   SetUpLCDVTopl(track_list,ip,bfld,pjet);
00057 }
00058 //__________________________________________________________________________
00059 
00060 //__________________________________________________________________________
00061 LCDVTopl::~LCDVTopl(){
00062   //  fclose(fp);
00063 
00064   if (vtrack_list) {
00065     vtrack_list->Delete();
00066     delete vtrack_list;
00067     vtrack_list=0;
00068   }
00069   if (vertex_list) {
00070     vertex_list->Delete();
00071     delete vertex_list;
00072     vertex_list=0;
00073   }
00074   if (vrttrack_list) {
00075     delete vrttrack_list;
00076     vrttrack_list=0;
00077   }
00078 }
00079 //__________________________________________________________________________
00080 
00081 //__________________________________________________________________________
00082 void LCDVTopl::SetUpLCDVTopl(TObjArray* track_list,
00083                              Int_t ntrk, 
00084                              Int_t* track_indecies,
00085                              const TVector3& ip,
00086                              const TVector3& bfld) {
00087   Clean();
00088 
00089   ipv   =ip;
00090   bfield=bfld;
00091   SetLCDVToplTRKs(track_list,ntrk,track_indecies);
00092   //This loop to 
00093   // make sure position and momentum in trk at 2d POCA
00094   // and  to make LCDVToplTRK.
00095   ptks.SetXYZ(0.0,0.0,0.0); 
00096   Double_t     l  =0.0;
00097   LCDVToplTRK* trk;
00098   for (Int_t i=0; i < ntrk ; i++) {
00099     trk= GetTrackAt(i);
00100     //Accumurate track momentum.
00101     ptks += trk->GetMomentumVector(l);
00102   } //loop over zxtrks
00103   pja=ptks;
00104 
00105   CalcIPErrorMatrix();
00106 }
00107 //__________________________________________________________________________
00108 
00109 //__________________________________________________________________________
00110 void LCDVTopl::SetUpLCDVTopl(TObjArray* track_list,
00111                              Int_t ntrk, 
00112                              Int_t* track_indecies,
00113                              const TVector3& ip,
00114                              const TVector3& bfld,
00115                              const TVector3& pjet) {
00116   Clean();
00117 
00118   ipv   =ip;
00119   bfield=bfld;
00120   SetLCDVToplTRKs(track_list,ntrk,track_indecies);
00121   //This loop to 
00122   // make sure position and momentum in trk at 2d POCA
00123   // and  to make LCDVToplTRK.
00124   ptks.SetXYZ(0.0,0.0,0.0); 
00125   Double_t     l  =0.0;
00126   LCDVToplTRK* trk;
00127   for( Int_t i=0; i < ntrk ; i++ ) {
00128     trk= GetTrackAt(i);
00129     //Accumurate track momentum.
00130     ptks += trk->GetMomentumVector(l);
00131   } //loop over zxtrks
00132   pja=pjet;
00133 
00134   CalcIPErrorMatrix();
00135 }
00136 //__________________________________________________________________________
00137 
00138 //__________________________________________________________________________
00139 void LCDVTopl::SetUpLCDVTopl(TObjArray* track_list,
00140                              const TVector3& ip,
00141                              const TVector3& bfld,
00142                              const TVector3& pjet) {
00143 
00144   Clean();
00145 
00146   ipv   =ip;
00147   bfield=bfld;
00148   LCDVToplTRK* trk;
00149   Int_t i;
00150   for(i=0 ; i < track_list->GetEntries() ; i++) {
00151     trk= new((*vtrack_list)[i]) LCDVToplTRK(track_list,i,ipv);
00152   } //loop over zxtrks
00153   //This loop to 
00154   // make sure position and momentum in trk at 2d POCA
00155   // and  to make LCDVToplTRK.
00156   ptks.SetXYZ(0.0,0.0,0.0); 
00157   Double_t     l  =0.0;
00158   for(i=0; i < track_list->GetEntries() ; i++ ) {
00159     trk= GetTrackAt(i);
00160     //Accumurate track momentum.
00161     ptks += trk->GetMomentumVector(l);
00162   } //loop over zxtrks
00163   pja=pjet;
00164 
00165   CalcIPErrorMatrix();
00166 }
00167 //__________________________________________________________________________
00168 
00169 //__________________________________________________________________________
00170 //Init...
00171 void LCDVTopl::Init() {
00172   iper.ResizeTo(3,3);
00173   iper_inv.ResizeTo(3,3);
00174   vertex_list   = new TClonesArray("LCDVToplVRT",1000);
00175   vtrack_list   = new TClonesArray("LCDVToplTRK",1000);
00176   vrttrack_list = new TObjArray(100);
00177     
00178   InitControlParameter();
00179 
00180   dblv=0;
00181 
00182   // Initialize jet total momentum
00183   pja.SetXYZ(0.0,0.0,0.0); 
00184 
00185   ptks.SetXYZ(0.0,0.0,0.0); 
00186 
00187   //Set Beam Position (Temporary)  
00188   ipv.SetXYZ(0.0,0.0,0.0); 
00189 
00190 }
00191 //__________________________________________________________________________
00192 
00193 //__________________________________________________________________________
00194 //Delete LCDVToplVRT, LCDVToplTRK...
00195 void LCDVTopl::Clean() {
00196   if (vtrack_list)   vtrack_list->Delete();
00197   if (vertex_list)   vertex_list->Delete();
00198   if (vrttrack_list) vrttrack_list->Clear();
00199 }
00200 //__________________________________________________________________________
00201 
00202 //__________________________________________________________________________
00203 //Get IP probabirity
00204 Double_t LCDVTopl::GetIPProbabirity(TVector3 xvtx) {
00205   TMatrixD ip2xvtx(3,1);
00206   for (Int_t i=0 ; i < 3 ; i++) { ip2xvtx(i,0)=xvtx[i] - ipv[i]; }
00207   TMatrixD tmp0(iper_inv,TMatrixD::kMult,ip2xvtx);
00208   TMatrixD tmp1(ip2xvtx,TMatrixD::kTransposeMult,tmp0);
00209   Double_t prfu = 0.0;
00210   // protect against floating overflow
00211   if (tmp1(0,0) < 100.0) {
00212     prfu = TMath::Exp(-0.5*tmp1(0,0));
00213   }
00214 
00215   return prfu;
00216 }
00217 //__________________________________________________________________________
00218 
00219 //__________________________________________________________________________
00220 //Get Initial vertexposition rrack and track 
00221 TVector3 LCDVTopl::GetInitialVertexT2T(Int_t itrk1, Int_t itrk2) {
00222   //initial vertex point for given tracks
00223 
00224   LCDVToplTRK* trk1= GetTrackAt(itrk1);
00225   TVector3 posv1   = trk1->GetPositionVector(0.0);
00226   // Double_t phi0pr1 = trk1->GetTrackParameter(1)-TMath::Pi();
00227   Double_t phi0pr1 = trk1->GetTrackParameter(1)-TMath::Pi()/2.0;
00228   Double_t cspr1   = TMath::Cos(phi0pr1);
00229   Double_t snpr1   = TMath::Sin(phi0pr1);
00230   Double_t z0pr1   = trk1->GetTrackParameter(3);
00231   Double_t x0pr1   = cspr1*(posv1[0] - ipv[0]) + snpr1*(posv1[1] - ipv[1]);
00232   Double_t tanlpr1 = trk1->GetTrackParameter(4);
00233   Double_t s1spr1  = trk1->GetErrorTsiApprox(0.0);
00234   Double_t s2spr1  = trk1->GetErrorEtaApprox(0.0)*(1.0 + tanlpr1*tanlpr1);
00235 
00236   LCDVToplTRK* trk2= GetTrackAt(itrk2);
00237   TVector3 posv2   = trk2->GetPositionVector(0.0);
00238   // Double_t phi0pr2 = trk2->GetTrackParameter(1)-TMath::Pi();
00239   Double_t phi0pr2 = trk2->GetTrackParameter(1)-TMath::Pi()/2.0;
00240   Double_t cspr2   = TMath::Cos(phi0pr2);
00241   Double_t snpr2   = TMath::Sin(phi0pr2);
00242   Double_t z0pr2   = trk2->GetTrackParameter(3);
00243   Double_t x0pr2   = cspr2*(posv2[0] - ipv[0]) + snpr2*(posv2[1] - ipv[1]);
00244   Double_t tanlpr2 = trk2->GetTrackParameter(4);
00245   Double_t s1spr2  = trk2->GetErrorTsiApprox(0.0);
00246   Double_t s2spr2  = trk2->GetErrorEtaApprox(0.0)*(1.0 + tanlpr2*tanlpr2);
00247 
00248   TMatrixD rot(3,3);
00249   rot(0,0) = cspr1*cspr1/s1spr1 + tanlpr1*tanlpr1*snpr1*snpr1/s2spr1 +
00250     cspr2*cspr2/s1spr2 + tanlpr2*tanlpr2*snpr2*snpr2/s2spr2;
00251   rot(1,0) = cspr1*snpr1/s1spr1 - tanlpr1*tanlpr1*snpr1*cspr1/s2spr1 +
00252     cspr2*snpr2/s1spr2 - tanlpr2*tanlpr2*snpr2*cspr2/s2spr2;
00253   rot(2,0) = tanlpr1*snpr1/s2spr1 + tanlpr2*snpr2/s2spr2;
00254   rot(0,1) = rot(1,0);
00255   rot(1,1) = snpr1*snpr1/s1spr1 + tanlpr1*tanlpr1*cspr1*cspr1/s2spr1 +
00256     snpr2*snpr2/s1spr2 + tanlpr2*tanlpr2*cspr2*cspr2/s2spr2;
00257   rot(2,1) = -(tanlpr1*cspr1/s2spr1 + tanlpr2*cspr2/s2spr2);
00258   rot(0,2) = rot(2,0);
00259   rot(1,2) = rot(2,1);
00260   rot(2,2) = 1/s2spr1 + 1/s2spr2;
00261   TMatrixD tmp(3,1);
00262   tmp(0,0) = -(cspr1*x0pr1/s1spr1 + z0pr1*tanlpr1*snpr1/s2spr1 +
00263                  cspr2*x0pr2/s1spr2 + z0pr2*tanlpr2*snpr2/s2spr2);
00264   tmp(1,0) = -snpr1*x0pr1/s1spr1 + z0pr1*tanlpr1*cspr1/s2spr1 -
00265     snpr2*x0pr2/s1spr2 + z0pr2*tanlpr2*cspr2/s2spr2;
00266   tmp(2,0) = -(z0pr1/s2spr1 + z0pr2/s2spr2);
00267   
00268   TMatrixD rot_inv(TMatrixD::kInverted,rot);
00269   TVector3 xdv=ipv;
00270   TMatrixD tmp1(3,1);
00271   tmp1.Mult(rot_inv,tmp);
00272   // xdv.SetXYZ(tmp1(0,0)+ipv[0],tmp1(1,0)+ipv[1],tmp1(2,0));
00273   xdv.SetXYZ(-tmp1(0,0)+ipv[0],-tmp1(1,0)+ipv[1],-tmp1(2,0));
00274   // avoid being too far off
00275   TVector3 xdv2(xdv[0],xdv[1],xdv[2]-ipv[2]);
00276   if ( xdv2.Mag2() > 16.0 ) {
00277       xdv=ipv;
00278   }
00279 
00280   return xdv;
00281 }
00282 //__________________________________________________________________________
00283 
00284 //__________________________________________________________________________
00285 //Get Initial vertexposition ip and track 
00286 TVector3 LCDVTopl::GetInitialVertexIP2T(Int_t itrk) {
00287   //initial vertex point for given tracks
00288   LCDVToplTRK* trk=GetTrackAt(itrk);
00289   TVector3 xpt=trk->GetPositionVector(0.0);
00290   Double_t stk2=trk->GetErrorDistance(0.0,ipv[0],ipv[1],ipv[2]);
00291   TVector3 dx=xpt-ipv;
00292   Double_t adx=dx.Mag();
00293   TMatrixD dddxT(3,1);
00294   for (Int_t i=0 ; i < 3 ; i++) { dddxT(i,0)=dx[i]/adx; }
00295   TMatrixD tmp0(iper,TMatrixD::kMult,dddxT);
00296   TMatrixD tmp1(dddxT,TMatrixD::kTransposeMult,tmp0);
00297   Double_t sip2=tmp1(0,0);
00298   Double_t wtot=1.0/sip2 + 1.0/stk2;
00299   return ipv + (1.0/wtot/stk2)*dx;
00300 
00301 }
00302 //__________________________________________________________________________
00303 
00304 //__________________________________________________________________________
00305 //Get vertex significance
00306 Double_t LCDVTopl::GetVertexSignificance(TVector3 xvrt) {
00307   //return vertex significance V(r) for NTRK tracks at x,y,z
00308  
00309   TVector3 ip2xvrt=xvrt-ipv;
00310 
00311   Double_t dlong=pja.Mag()> 0 ? (ip2xvrt*pja)/pja.Mag() : ip2xvrt.Mag();
00312 
00313   // cut out if behind ip or > 10.0 cm in front
00314   if (dlong > 10.0 || dlong < -0.01) {
00315     return -1.0;
00316   }
00317 
00318   Double_t dmag = ip2xvrt.Mag();
00319   Double_t dtran= 0.0;
00320   if (TMath::Abs(dlong) < TMath::Abs(dmag)) {
00321     dtran = TMath::Sqrt(dmag*dmag - dlong*dlong); 
00322   }
00323  
00324   Double_t prsumu  = 0.0;
00325   Double_t prsum2u = 0.0;
00326   Double_t prfu=0.0;
00327   Int_t ntrk=GetTrackEntries();
00328   LCDVToplTRK* trk=0;
00329   Int_t i=0;
00330   for (i=0 ; i < ntrk ; i++) { 
00331     trk =GetTrackAt(i);
00332     if (trk->GetTrackID() < 0) { continue; } //skip IP track...
00333     //prfu=trk->GetTrackProbabirity(xvrt);
00334     prfu=trk->GetTrackProbabirity(xvrt,ipv);
00335     prsumu  += prfu;
00336     prsum2u += prfu*prfu;
00337   } //End of the track loop.
00338   
00339   // add in IP information if required
00340   prfu = GetIPProbabirity(xvrt);
00341   
00342   prsumu  += kipw*prfu;
00343   prsum2u += kipw*prfu*prfu; 
00344         
00345   Double_t vsg = 0.0;
00346   if (prsumu > 0.0) {
00347     vsg = prsumu - prsum2u/prsumu; 
00348   };
00349  
00350   //include extra factors such as jet core weighting
00351   // What are 0.01 and 0.005?
00352   if (dtran > 0.005) {
00353     Double_t alpha = TMath::ACos((dlong + 0.01)/
00354                                  TMath::Sqrt(TMath::Power(dlong + 0.01,2)
00355                                              + TMath::Power(dtran - 0.005,2)));
00356     //What is 0.01 and 0.005????
00357     vsg *= TMath::Exp(-kang*alpha*alpha);
00358   }
00359 
00360   return TMath::Max(vsg,0.0);
00361 }
00362 //__________________________________________________________________________
00363 
00364 //__________________________________________________________________________
00365 //Get min(Vr/min(V))
00366 Double_t LCDVTopl::GetMinVrOverMinVs(LCDVToplVRT* vrt1, LCDVToplVRT* vrt2, Int_t nstp) {
00367 
00368   Double_t v1   =vrt1->GetVsig();
00369   Double_t v2   =vrt2->GetVsig();
00370   Double_t vmin =TMath::Min(v1,v2);
00371   TVector3 xvrt1=vrt1->GetVertexVector();
00372   TVector3 xvrt2=vrt2->GetVertexVector();
00373   TVector3 dxvrt=xvrt2-xvrt1;
00374 
00375   // prevent having to find vmin if vertices are within 10 um
00376   Double_t dstep =dxvrt.Mag();
00377   Double_t vrat  =0;
00378   Double_t vmaxi =0;
00379   Double_t vmini =0;
00380   Int_t    mstep =nstp;
00381   TVector3 xpt(0.0,0.0,0.0);
00382   if (dstep < dcut) {
00383     vrat = 1.0;
00384   } else {
00385     dstep /= nstp + 1;
00386     if (dstep < dcut/2.0) {
00387       mstep=(Int_t)((nstp+1)*dstep/dcut*2.0) - 1;
00388       dstep = (nstp + 1)*dstep/(mstep + 1);
00389     }
00390     dxvrt.SetMag(dstep);
00391     vmini=vmin;
00392     mstep=TMath::Max(mstep,1);
00393     for (Int_t step=1 ; step <= mstep ; step++) {
00394       xpt  =xvrt1 + step*dxvrt;
00395       vmaxi=GetVertexSignificance(xpt);
00396       if (vmaxi < vmini) { vmini = vmaxi; }
00397     } //endloop over steps
00398     vrat = vmini/vmin;
00399   }
00400 
00401   return vrat;
00402 }
00403 //__________________________________________________________________________
00404 
00405 //__________________________________________________________________________
00406 Double_t LCDVTopl::GetMinVrOverMinVs(LCDVToplVRT* vrt1, LCDVToplVRT* vrt2) {
00407   return GetMinVrOverMinVs(vrt1,vrt2,nstep);
00408 }
00409 //__________________________________________________________________________
00410 
00411 //__________________________________________________________________________
00412 //Get Vertex by simple iteration
00413 TVector3 LCDVTopl::GetVertexBySimpleIteration(LCDVToplVRT* vrt) {
00414 
00415   // ************ 3D iterative bit ****************
00416 
00417   TVector3 spx       = vrt->GetVertexVector();
00418   Double_t vm3d      = vrt->GetVsig();
00419   //start with 4 um step size
00420   Double_t stepsize     = 0.0005;
00421   Bool_t f_iteration = kFALSE;
00422   Bool_t f_again     = kFALSE;
00423   Double_t vsg       = 0.0;
00424 
00425   while(kTRUE) {
00426 
00427     f_iteration = kFALSE;
00428 
00429     while(1) {
00430       spx += TVector3(stepsize,0.0,0.0);
00431       vsg=GetVertexSignificance(spx);
00432       if (vsg > vm3d) {
00433         vm3d = vsg;
00434         f_iteration = kTRUE;
00435       } else {
00436         spx -= TVector3(stepsize,0.0,0.0);
00437         break;
00438       }
00439     }
00440  
00441     if (!f_iteration) {
00442       while(1) {
00443         spx -= TVector3(stepsize,0.0,0.0);
00444         vsg=GetVertexSignificance(spx);
00445         if (vsg > vm3d) {
00446           vm3d = vsg;
00447           f_iteration = kTRUE;
00448         } else {
00449           spx += TVector3(stepsize,0.0,0.0);
00450           break;
00451         }
00452       }
00453     }
00454 
00455     f_again=kFALSE;
00456 
00457     //now y co-ordinate
00458  
00459     f_iteration = kFALSE;
00460 
00461     while(1) {
00462       spx += TVector3(0.0,stepsize,0.0);
00463       vsg=GetVertexSignificance(spx);
00464       if (vsg > vm3d) {
00465         f_again = kTRUE;
00466         vm3d  = vsg;
00467         f_iteration  = kTRUE;
00468       } else {
00469         spx -= TVector3(0.0,stepsize,0.0);
00470         break;
00471       }
00472     }
00473     if (!f_iteration) {
00474       while(1) {
00475         spx -= TVector3(0.0,stepsize,0.0);
00476         vsg=GetVertexSignificance(spx);
00477         if (vsg > vm3d) {
00478           f_again = kTRUE;
00479           vm3d  = vsg;
00480         } else {
00481           spx += TVector3(0.0,stepsize,0.0);
00482           break;
00483         }
00484       }
00485     }
00486     
00487     //now z co-ordinate
00488     
00489     f_iteration = kFALSE;
00490     while(1) {
00491       spx += TVector3(0.0,0.0,stepsize);
00492       vsg=GetVertexSignificance(spx);
00493       if (vsg > vm3d) {
00494         f_again = kTRUE;
00495         vm3d  = vsg; 
00496         f_iteration  = kTRUE;
00497       } else {
00498         spx -= TVector3(0.0,0.0,stepsize);
00499         break;
00500       }
00501     }
00502     
00503     if (!f_iteration) {
00504       while(1) {
00505         spx -= TVector3(0.0,0.0,stepsize);
00506         vsg=GetVertexSignificance(spx);
00507         if (vsg > vm3d) {
00508           f_again = kTRUE;
00509           vm3d  = vsg; 
00510         } else {
00511           spx += TVector3(0.0,0.0,stepsize);
00512           break;
00513         }
00514       }
00515     }
00516     
00517     if (!f_again) {
00518       if (stepsize > dcut) {
00519         stepsize = stepsize/2.0;
00520       } else {
00521         break;
00522       }
00523     }
00524   }
00525 
00526   vrt->SetVsig(vm3d);
00527   vrt->SetVertex(spx);
00528   return spx;
00529 }
00530 //__________________________________________________________________________
00531 
00532 //__________________________________________________________________________
00533 //Get Vertex Charge
00534 Double_t LCDVTopl::GetVertexCharge() {
00535   Int_t nvtrks=GetVrtTrackEntries();
00536   LCDVToplTRK* trk=0;
00537   Double_t q_vertex=0;
00538   for (Int_t itrk=0 ; itrk < nvtrks ; itrk++) {
00539     trk=GetVrtTrackAt(itrk);
00540     q_vertex += trk->GetCharge();
00541   }
00542   return q_vertex;
00543 }
00544 //__________________________________________________________________________
00545 
00546 //__________________________________________________________________________
00547 //Get # of tracks
00548 Int_t LCDVTopl::GetTrackEntries(){
00549     return vtrack_list->GetEntries();
00550 }
00551 //__________________________________________________________________________
00552 
00553 //__________________________________________________________________________
00554 LCDVToplTRK* LCDVTopl::GetTrackAt(Int_t itrk) {
00555   return (LCDVToplTRK*)(vtrack_list->At(itrk));
00556 }
00557 //__________________________________________________________________________
00558 
00559 //__________________________________________________________________________
00560 Int_t LCDVTopl::GetVertexEntries(){
00561     return vertex_list->GetEntries();
00562 }
00563 //__________________________________________________________________________
00564 
00565 //__________________________________________________________________________
00566 Int_t LCDVTopl::GetVertexEntriesExceptV0(){
00567   Int_t nvrt=GetVertexEntries();
00568   Int_t nvrtnov0=0;
00569   LCDVToplVRT* vrt;
00570   for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
00571     vrt=GetVertexAt(ivrt);
00572     if (CheckKs0Vertex(vrt))     continue;
00573     if (CheckLambda0Vertex(vrt)) continue;
00574     nvrtnov0++;
00575   }
00576   return nvrtnov0;
00577 }
00578 //__________________________________________________________________________
00579 
00580 //__________________________________________________________________________
00581 LCDVToplVRT* LCDVTopl::GetVertexAt(Int_t ivrt) {
00582   return (LCDVToplVRT*)(vertex_list->At(ivrt));
00583 }
00584 //__________________________________________________________________________
00585 
00586 //__________________________________________________________________________
00587 Int_t LCDVTopl::GetVrtTrackEntries(){
00588     return vrttrack_list->GetEntries();
00589 }
00590 //__________________________________________________________________________
00591 
00592 //__________________________________________________________________________
00593 LCDVToplTRK* LCDVTopl::GetVrtTrackAt(Int_t itrk) {
00594   return (LCDVToplTRK*)(vrttrack_list->At(itrk));
00595 }
00596 //__________________________________________________________________________
00597 
00598 //__________________________________________________________________________
00599 void LCDVTopl::InitControlParameter(){
00600   //   ripe  =   3.0;
00601   //   zipe  =  10.0;
00602    ripe  =  3.0e-4;
00603    zipe  =  6.0e-4;
00604    rcut  =  0.6;
00605    xcut  = 10.0;
00606    kang  =  5.0;
00607    kipw  =  1.0;
00608 
00609    nstep = 7;
00610    dcut  = 0.0002;
00611 
00612    m_SeedVrtRMax  =2.3  ;
00613    m_SeedVrtK0Mass=0.015;
00614    m_SeedVrtLambda0Mass=0.015;
00615    m_SeedVrtSigMin=2.0  ;
00616 
00617    m_TrdiMax=0.1;
00618    m_LodiMin=0.05;
00619    m_LodrMin=0.25;
00620    //m_LodrMax=2.0;
00621    m_LodrMax=100.0;
00622    //m_AntaMax=0.5;
00623    m_AntaMax=TMath::Pi()/2.0;
00624 
00625    //m_PTIP[0]=0.0010;
00626    //m_PTIP[1]=0.0040;
00627    m_PTIP[0]=1.5*ripe;
00628    m_PTIP[1]=1.5*zipe;
00629    m_NVsigma=1.0;
00630 
00631    m_Rmass   =0.0;
00632    m_PtMass  =0.0;
00633    m_PtMassEr=0.0;
00634    m_Masstag =0.0;
00635    m_MissPt  =0.0;
00636 
00637    m_UseNN     =kFALSE;
00638    m_SeedVrtNN =0.7;
00639    m_TrackAttNN=0.4;
00640    m_SelectNN  =0.5;
00641 
00642    m_p4Vertex.SetXYZT(0.,0.,0.,0.);
00643 }
00644 //__________________________________________________________________________
00645 
00646 //__________________________________________________________________________
00647 void LCDVTopl::SetLCDVToplTRKs(TObjArray* track_list,
00648                               Int_t ntrk, 
00649                               Int_t* track_indecies) {
00650   //This loop to 
00651   // make sure position and momentum in trk at 2d POCA
00652   // and  to make LCDVToplTRK.
00653   //  pja.SetXYZ(0.0,0.0,0.0); 
00654   for( Int_t i=0; i < ntrk ; i++ ) {
00655     new((*vtrack_list)[i]) LCDVToplTRK(track_list,track_indecies[i],ipv);
00656   } //loop over zxtrks
00657 }
00658 //__________________________________________________________________________
00659 
00660 //__________________________________________________________________________
00661 void LCDVTopl::CalcIPErrorMatrix(){
00662   //Double_t ripe2=ripe/1e8*ripe;
00663   //Double_t zipe2=zipe/1e8*zipe;
00664   Double_t ripe2=ripe*ripe;
00665   Double_t zipe2=zipe*zipe;
00666   iper(0,0)=iper_inv(0,0)=ripe2;
00667   iper(1,0)=iper_inv(1,0)=0.0  ;
00668   iper(2,0)=iper_inv(2,0)=0.0  ;
00669   iper(0,1)=iper_inv(0,1)=0.0  ;
00670   iper(1,1)=iper_inv(1,1)=ripe2;
00671   iper(2,1)=iper_inv(2,1)=0.0  ;
00672   iper(0,2)=iper_inv(0,2)=0.0  ;
00673   iper(1,2)=iper_inv(1,2)=0.0  ;
00674   iper(2,2)=iper_inv(2,2)=zipe2;
00675   Double_t det_iper;
00676   iper_inv.Invert(&det_iper);
00677 
00678 }
00679 //__________________________________________________________________________
00680 
00681 //__________________________________________________________________________
00682 void LCDVTopl::DeleteTracks(){
00683   vtrack_list->Delete();
00684 }
00685 //__________________________________________________________________________
00686 
00687 //__________________________________________________________________________
00688 void LCDVTopl::DeleteVertices(){
00689   vertex_list->Delete();
00690 }
00691 //__________________________________________________________________________
00692 
00693 //__________________________________________________________________________
00694 void LCDVTopl::AddVertex(LCDVToplVRT* vrt) {
00695   Int_t nvrt=GetVertexEntries();
00696   Int_t ivrt=0;
00697   Bool_t f_no_same_vertex=kTRUE;
00698   for (ivrt=0 ; ivrt < nvrt-1 ; ivrt++) {
00699     if (vrt == GetVertexAt(ivrt)) {
00700       f_no_same_vertex=kFALSE;
00701       break;
00702     }
00703   }
00704 
00705   if (f_no_same_vertex) {
00706     SortVertexByVsig();
00707   } else {
00708     vertex_list->Remove(vrt);
00709   }
00710 }
00711 //__________________________________________________________________________
00712 
00713 //__________________________________________________________________________
00714 void LCDVTopl::AddVertexIntoTrack(LCDVToplTRK* trk, LCDVToplVRT* vrt) {
00715   Int_t nvrt=trk->GetVertexEntries();
00716   Int_t ivrt=0;
00717   Bool_t f_no_same_vertex=kTRUE;
00718   LCDVToplVRT* tvrt=0;
00719   for (ivrt=0; ivrt < nvrt ; ivrt++) {
00720     tvrt=(LCDVToplVRT*)trk->GetVertexAt(ivrt);
00721     if (tvrt == vrt) { 
00722       f_no_same_vertex=kFALSE; 
00723       break;
00724     }
00725   }
00726   if (f_no_same_vertex) {
00727     trk->AddVertex(vrt);
00728     SortVertexByVsigInTrack(trk);
00729   }  
00730 }
00731 //__________________________________________________________________________
00732 
00733 //__________________________________________________________________________
00734 void LCDVTopl::AddVrtTrack(LCDVToplTRK* trk) {
00735   vrttrack_list->Add(trk);
00736 }
00737 //__________________________________________________________________________
00738 
00739 //__________________________________________________________________________
00740 void LCDVTopl::CompressVerticesList() {
00741   vertex_list->Compress();
00742 }
00743 //__________________________________________________________________________
00744 
00745 //__________________________________________________________________________
00746 void LCDVTopl::RemoveVertexAt(Int_t i) {
00747   RemoveVertex(GetVertexAt(i));
00748 }
00749 //__________________________________________________________________________
00750 
00751 //__________________________________________________________________________
00752 void LCDVTopl::RemoveVertex(LCDVToplVRT* vrt) {
00753   vertex_list->Remove(vrt);
00754   CompressVerticesList();
00755 }
00756 //__________________________________________________________________________
00757 
00758 //__________________________________________________________________________
00759 void LCDVTopl::MergeVertex(LCDVToplVRT* vrti,LCDVToplVRT* vrtj,Int_t f_fit) {
00760   if (vrti == vrtj) { return; }
00761 
00762   LCDVToplTRK* trk=0;
00763   Int_t ntrk=vrtj->GetTrackEntries();
00764   Int_t itrk=0;
00765   for (itrk=0; itrk < ntrk ; itrk++) {
00766     trk =vrtj->GetTrackAt(itrk);
00767     vrti->AddTrack(trk,f_fit);
00768     trk->RemoveVertex(vrtj);
00769     AddVertexIntoTrack(trk,vrti);
00770     trk->CompressVerticesList();
00771   }
00772 
00773   if (f_fit > 1) {
00774     vrti->SetVsig(GetVertexSignificance(vrti->GetVertexVector()));
00775   }
00776 
00777   RemoveVertex(vrtj);
00778 
00779   CompressVerticesList();
00780 }
00781 //__________________________________________________________________________
00782 
00783 //__________________________________________________________________________
00784 void LCDVTopl::MergeVertex(Int_t i, Int_t j, Int_t f_fit) {
00785   if (i == j) { return; }
00786 
00787   LCDVToplVRT* vrti=GetVertexAt(i);
00788   LCDVToplVRT* vrtj=GetVertexAt(j);
00789 
00790   MergeVertex(vrti,vrtj,f_fit);
00791 }
00792 //__________________________________________________________________________
00793 
00794 //__________________________________________________________________________
00795 void LCDVTopl::SortVertexByVsig() {
00796   Int_t nvrt = vertex_list->GetEntries();
00797   for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
00798     LCDVToplVRT* vrt=(LCDVToplVRT*)vertex_list->At(ivrt);
00799     Double_t vsig=vrt->GetVsig();
00800     vrt->SetSortVal(-vsig);
00801   }
00802   vertex_list->Sort();
00803 }
00804 //__________________________________________________________________________
00805 
00806 //__________________________________________________________________________
00807 void LCDVTopl::SortVertexByDecayLength() {
00808   Int_t nvrt=GetVertexEntries();
00809   Int_t    ivrti,ivrtj;
00810   Double_t disti,distj;
00811   LCDVToplVRT* vrti;
00812   LCDVToplVRT* vrtj;
00813   LCDVToplVRT  vrtk;
00814   for (ivrti=0 ; ivrti < nvrt-1 ; ivrti++) {
00815     for (ivrtj=ivrti+1 ; ivrtj < nvrt ; ivrtj++) {
00816       vrti=GetVertexAt(ivrti);
00817       vrtj=GetVertexAt(ivrtj);
00818       disti=vrti->GetDecayLength(ipv);
00819       distj=vrtj->GetDecayLength(ipv);
00820       if (disti > distj) {
00821         vrtk = *vrtj;
00822         *vrtj = *vrti;
00823         *vrti = vrtk;
00824       }
00825     }
00826   }
00827 }
00828 //__________________________________________________________________________
00829 
00830 //__________________________________________________________________________
00831 void LCDVTopl::SortVertexByVsigInTrack(LCDVToplTRK* trk) {
00832   Int_t nvrt=trk->GetVertexEntries();
00833   Int_t ivrti=0,ivrtj=0;
00834   LCDVToplVRT* vrti=0;
00835   LCDVToplVRT* vrtj=0;
00836 
00837   for (ivrti=0 ; ivrti < nvrt-1 ; ivrti++) {
00838     for (ivrtj=ivrti+1 ; ivrtj < nvrt ; ivrtj++) {
00839       vrti=(LCDVToplVRT*)(trk->GetVertexAt(ivrti));
00840       vrtj=(LCDVToplVRT*)(trk->GetVertexAt(ivrtj));
00841       if (vrti->GetVsig() < vrtj->GetVsig()) {
00842         // trk->AddVertexAt((TObject*)vrtj,ivrti);
00843         // trk->AddVertexAt((TObject*)vrti,ivrtj);
00844         trk->AddVertexAt(vrtj,ivrti);
00845         trk->AddVertexAt(vrti,ivrtj);
00846       }
00847     }
00848   }
00849 
00850 }
00851 //__________________________________________________________________________
00852 
00853 //__________________________________________________________________________
00854 void LCDVTopl::SortVertexByVsigInAllTrack() {
00855   Int_t ntrk=GetTrackEntries();
00856   for (Int_t itrk=0 ; itrk < ntrk ; itrk++) {
00857     SortVertexByVsigInTrack(GetTrackAt(itrk));
00858   }
00859 }
00860 //__________________________________________________________________________
00861 
00862 //__________________________________________________________________________
00863 Int_t LCDVTopl::FindVertices() {
00864   return FindVertices(0);
00865 }
00866 //__________________________________________________________________________
00867 
00868 //__________________________________________________________________________
00869 Int_t LCDVTopl::FindVertices(Int_t f_v0check) {
00870   
00871   //-------------------------------------------------------
00872   //
00873   // FindVertices --  
00874   //    determine the TOPOlogical Vertex structure of a jet
00875   //------------------------------------------------------
00876   //
00877   // resolve_algo -- resolve track V(r) into inclusive vertices
00878   //-----------------------------------------------------
00879 
00880   if (dblv > 0)
00881     cout << "=========<FindVertices started!>============"<<endl;
00882   
00883   Int_t ntrk=GetTrackEntries();
00884   
00885   if (dblv > 0)
00886     cout << "   ---> Check 1: N of LCDVToplTRKs= " << ntrk << endl;
00887 
00888   // first loop over track pairs to find Vs and interactions
00889   if (ntrk < 2) {
00890     //not enough tracks!
00891     return 0;
00892   }
00893   
00894   if (dblv > 0)
00895     cout <<"   ---> Check 2: CalcIPErrorMatrix() ... ";
00896   CalcIPErrorMatrix();
00897   if (dblv > 0)
00898     cout <<"done"<<endl;
00899 
00900   // Addd IP-Track...
00901   if (dblv > 0)
00902     cout << "   ---> Check 3: AddIPTrack ... ";
00903   LCDVToplTRK* iptrk=new((*vtrack_list)[ntrk]) LCDVToplTRK();
00904   if (dblv > 0)
00905     cout <<" done"<<endl;
00906 
00907   //----------------------------------------------------
00908   // now loop over track pairs and IP
00909   //----------------------------------------------------
00910   if (dblv > 0)
00911     cout << "   ---> Check 4: Making Vertices with IP-Track ... ";
00912   TVector3 xvinit(0.0,0.0,0.0);
00913   TVector3 xvrt(0.0,0.0,0.0);
00914   LCDVToplTRK* trkj = 0;
00915   LCDVToplTRK* trk  = 0;
00916   LCDVToplVRT* vrt  = 0;
00917   Double_t     vmaxi= 0.0;
00918   Double_t     prfi = 0;
00919   Double_t     prfj = 0;
00920   Double_t     lprfi= 0;
00921   Double_t     lprfj= 0;
00922   Int_t        itrk = ntrk;
00923   Int_t        jtrk = 0;
00924   Int_t        nvrt = 0;
00925   for (jtrk=0 ; jtrk < ntrk ;jtrk++) { 
00926     trkj=GetTrackAt(jtrk);
00927 
00928     //ip-track combination
00929     xvinit=GetInitialVertexIP2T(jtrk);
00930 
00931     vrt=new((*vertex_list)[nvrt]) LCDVToplVRT(xvinit);
00932     vrt->AddTrack(iptrk);
00933     vrt->AddTrack(trkj);
00934     vrt->SetVertex(xvinit);
00935 
00936     xvrt=vrt->GetVertexVector();
00937 
00938     vmaxi=GetVertexSignificance(xvrt);
00939 
00940     //Vertex significance check.
00941     if (vmaxi <= 0.001) {
00942       vertex_list->Remove(vrt);
00943       continue;
00944     }
00945 
00946     vrt->SetVsig(vmaxi);
00947     prfi=GetIPProbabirity(xvrt);
00948     prfj=trkj->GetTrackProbabirity(xvrt,ipv);
00949 
00950     //Track Probabirity Check.
00951     if (prfi*prfj <= 0.0) {
00952       vertex_list->Remove(vrt);
00953       continue;
00954     }
00955 
00956     lprfi = -2.0*TMath::Log(prfi);
00957     lprfj = -2.0*TMath::Log(prfj);
00958 
00959     //Chi**2 contribution of tracks check.
00960     if (lprfi > xcut || lprfj > xcut) {
00961       vertex_list->Remove(vrt);
00962       continue;
00963     }
00964 
00965     nvrt++;
00966     AddVertexIntoTrack(iptrk,vrt);
00967     AddVertexIntoTrack(trkj ,vrt);
00968 
00969   } // end do loop over track pairs
00970  
00971   //Int_t nvrt=GetVertexEntries();
00972   if (dblv > 0)
00973     cout <<nvrt<<" done"<<endl;
00974 
00975   //----------------------------------------------------
00976   // now loop over track pairs
00977   //----------------------------------------------------
00978   if (dblv > 0)
00979     cout << "   ---> Check 5: Making Vertices by Track-Track ... ";
00980   LCDVToplTRK* trki;
00981   for (itrk=1 ; itrk < ntrk ; itrk++) { 
00982     trki=GetTrackAt(itrk);
00983 
00984     for (jtrk=0 ; jtrk < itrk ;jtrk++) { 
00985 
00986       trkj=GetTrackAt(jtrk);
00987       
00988       //track-track combination
00989       xvinit=GetInitialVertexT2T(itrk,jtrk);
00990 
00991       vrt=new((*vertex_list)[nvrt]) LCDVToplVRT(xvinit);
00992       vrt->AddTrack(trki);
00993       vrt->AddTrack(trkj);
00994 
00995       xvrt=vrt->GetVertexVector();
00996 
00997       prfi=trki->GetTrackProbabirity(xvrt,ipv);
00998       prfj=trkj->GetTrackProbabirity(xvrt,ipv);
00999 
01000       //Track Probabirity Check.
01001       if (prfi*prfj <= 0.0) {
01002         vertex_list->Remove(vrt);
01003         continue;
01004       }
01005 
01006       lprfi = -2.0*TMath::Log(prfi);
01007       lprfj = -2.0*TMath::Log(prfj);
01008 
01009       //Chi**2 contribution of tracks check.
01010       if (lprfi > xcut || lprfj > xcut ) {
01011         vertex_list->Remove(vrt);
01012         continue;
01013       }
01014 
01015       vmaxi=GetVertexSignificance(xvrt);
01016       vrt->SetVsig(vmaxi);
01017 
01018       //Vertex significance check.
01019       if (vmaxi <= 0.001) {
01020         vertex_list->Remove(vrt);
01021         continue;
01022       }
01023 
01024       //V0 check...
01025       if (f_v0check != 0 && vrt->GetChi2() < 16.0) {
01026         vrt->CalcFourVector();
01027         if (CheckKs0Vertex(vrt)) {
01028           trki->SetV0flag(1);
01029           trkj->SetV0flag(1);
01030         }
01031         if (CheckLambda0Vertex(vrt)) {
01032           trki->SetV0flag(2);
01033           trkj->SetV0flag(2);
01034         }
01035       }
01036 
01037       //            xvrt=GetVertexBySimpleIteration(vrt);
01038       nvrt++;
01039       AddVertexIntoTrack(trki,vrt);
01040       AddVertexIntoTrack(trkj,vrt);
01041 
01042     } // end do loop over track pairs
01043 
01044   } //end do loop over track pairs
01045  
01046   SortVertexByVsig();
01047   SortVertexByVsigInAllTrack();
01048   //nvrt=GetVertexEntries();
01049   if (dblv > 0)
01050     cout <<nvrt<<" done"<<endl;
01051 
01052   Int_t ivrt;
01053   Int_t nvtrk = 0;
01054   PrintOutForDebug();
01055 
01056   // -----------> now resolve vmax for each track
01057   if (dblv > 0)
01058     cout << "   ---> Check 6: Resolve Vertices per track ... ";
01059   Double_t gmax= 0.0;
01060   Double_t vmax=-1.0;
01061   Int_t    itsp=0,ntsp=0;
01062   Int_t    irsp=0,nrsp=0;
01063   TVector3 xvrt_prev(0.0,0.0,0.0);
01064   TVector3 dvrt(0.0,0.0,0.0);
01065   Bool_t   f_resolved=kTRUE;
01066   LCDVToplVRT* vrt_prev=0;
01067   Double_t vrat=0.0;
01068 
01069   for (itrk=0 ; itrk < ntrk+1 ; itrk++ ) {
01070     //loop over tracks + ip
01071     trk=GetTrackAt(itrk);
01072     gmax = 0.0; // global max on track
01073     vmax =-1.0;
01074     nrsp = 0;
01075     ntsp = trk->GetVertexEntries();
01076 
01077     itsp=0;
01078     while (itsp < ntsp) { //loop over track's spatial points
01079       vrt=(LCDVToplVRT*)(trk->GetVertexAt(itsp));
01080       vmax = vrt->GetVsig();
01081 
01082       if (nrsp == 0) { gmax = vmax; }
01083 
01084       if (vmax > 0.001 && vmax > 0.1*gmax) { // must be >frac of big max
01085 
01086         xvrt=vrt->GetVertexVector();
01087         f_resolved=kTRUE;
01088         // no. of res. spatial points on track
01089         for (irsp=0 ; irsp < nrsp ; irsp++) {     
01090           vrt_prev=(LCDVToplVRT*)(trk->GetVertexAt(irsp));
01091           xvrt_prev=vrt_prev->GetVertexVector();
01092           dvrt=xvrt - xvrt_prev;
01093           if (dvrt.Mag() < dcut) {
01094             f_resolved=kFALSE;
01095             break;
01096           }
01097 
01098           vrat    =GetMinVrOverMinVs(vrt,vrt_prev,7);
01099 //        vrat    =GetMinVrOverMinVs(vrt,vrt_prev,15);
01100           if (vrat > rcut) {
01101             f_resolved=kFALSE; 
01102             break;
01103           }
01104 
01105         } //endloop over resolved max.
01106         if (f_resolved) {
01107           nrsp++;
01108           itsp++;
01109         } else {
01110           MergeVertex(vrt_prev,vrt,0);
01111           ntsp=trk->GetVertexEntries();
01112         } //endif f_resolved
01113 
01114       } else {
01115         trk->RemoveVertexAt(itsp);
01116         vrt->RemoveTrack(trk,0);
01117         ntsp=trk->GetVertexEntries();
01118       } // end if vmax > 0.001 && vmax > 0.1*gmax
01119 
01120     } //endloop over spatial points on track.
01121 
01122   } //endloop over tracks
01123 
01124   nvrt=GetVertexEntries();
01125   if (dblv > 0)
01126     cout <<nvrt<<" done"<<endl;
01127   PrintOutForDebug();
01128 
01129   //fill resolved spatial max. arrays
01130   if (dblv > 0)
01131     cout << "   ---> Check 7: Resolve Vertices by Track info. ";
01132   Int_t        nsp   = GetVertexEntries();
01133   Int_t        isp   = 0;
01134   Int_t        nresp = 0;
01135   Int_t        iresp = 0;
01136   LCDVToplVRT* tvrt  = 0;
01137   Bool_t       f_not_found_yet=kTRUE;
01138   //  Int_t k=0;
01139   isp=0;
01140   while(isp < nsp) {
01141     //  for (isp=0 ; isp < nsp ; isp++) {
01142     vrt=GetVertexAt(isp);
01143     nvtrk=vrt->GetTrackEntries();
01144     if (nvtrk > 1) {
01145       
01146       //use logical f_resolved for convenience
01147       f_resolved=kFALSE; 
01148       
01149       for (itrk=0 ; itrk < ntrk+1 ; itrk++) {
01150         trk =GetTrackAt(itrk);
01151         nrsp=trk->GetVertexEntries();
01152         for (irsp=0 ; irsp < nrsp ; irsp++) {
01153           tvrt=(LCDVToplVRT*)(trk->GetVertexAt(irsp));
01154           if (vrt == tvrt) {
01155             f_resolved=kTRUE;
01156             break;
01157           }
01158         }
01159         if (f_resolved) { break; }
01160       }
01161 
01162       if (f_resolved &&  vrt->GetVsig() > 0.001) {
01163         
01164         xvrt=GetVertexBySimpleIteration(vrt);
01165 
01166         //has this spatial max already been found
01167         f_not_found_yet=kTRUE;
01168         for (iresp=0 ; iresp < nresp ; iresp++) {
01169           vrt_prev =GetVertexAt(iresp);
01170           xvrt_prev=vrt_prev->GetVertexVector();
01171           dvrt     =xvrt - xvrt_prev;
01172           
01173           if (dvrt.Mag() < dcut) {
01174             //if old spatial point, assign tracks to it
01175             nvtrk=vrt->GetTrackEntries();
01176             for (itrk=0 ; itrk < nvtrk ; itrk++) { 
01177               trk =vrt->GetTrackAt(itrk);
01178               nrsp=trk->GetVertexEntries();
01179               for (irsp=0; irsp < nrsp ; irsp++) {
01180                 tvrt = (LCDVToplVRT*)(trk->GetVertexAt(irsp));
01181                 if (vrt == tvrt) {
01182                   trk->RemoveVertexAt(irsp);
01183                   AddVertexIntoTrack(trk,vrt_prev);
01184                   //T.Abe Add 3/16/2001
01185                   vrt_prev->AddTrack(trk);
01186                 }
01187               }
01188             }
01189             RemoveVertexAt(isp);
01190             nsp=GetVertexEntries();
01191             f_not_found_yet=kFALSE;
01192             break;
01193           }
01194 
01195         }
01196         
01197         if (f_not_found_yet) {
01198           nresp++;
01199           isp++;
01200         }
01201         
01202       } else {
01203         RemoveVertexAt(isp);
01204         nsp=GetVertexEntries();
01205       }// endif ....
01206 
01207     } else { // # of tracks in the vertex < 2
01208       if (nvtrk > 0) {
01209         trk=vrt->GetTrackAt(0);
01210         // trk->RemoveVertex((TObject*)vrt);
01211         trk->RemoveVertex(vrt);
01212       }
01213       RemoveVertexAt(isp);
01214       nsp=GetVertexEntries();
01215     } // endif nvtrk > 1
01216 
01217   } // endloop of isp
01218 
01219   nvrt=GetVertexEntries();
01220 
01221   SortVertexByVsig();
01222   SortVertexByVsigInAllTrack();
01223 
01224   if (dblv > 0)
01225     cout <<nvrt<<" done"<<endl;
01226   PrintOutForDebug();
01227 
01228   //
01229   // loop over spatial points, fill resolution matrix Vminsp
01230   //
01231 
01232   // now find resolved vertices within rcut
01233   if (dblv > 0)
01234     cout << "   ---> Check 8: Resolve Vertices by Vsig... ";
01235   nresp=GetVertexEntries();
01236   LCDVToplVRT* vrti=0;
01237   LCDVToplVRT* vrtj=0;
01238   iresp=0;
01239   Int_t jresp=0;
01240   while(iresp < nresp-1) {
01241     vrti=GetVertexAt(iresp);
01242     jresp = iresp + 1;
01243     while(jresp < nresp) {
01244       vrtj=GetVertexAt(jresp);
01245       vrat=GetMinVrOverMinVs(vrti,vrtj,7);
01246 //      vrat=GetMinVrOverMinVs(vrti,vrtj,15);
01247       if (vrat > rcut) {
01248         MergeVertex(iresp,jresp,0);
01249         nresp=GetVertexEntries();
01250       } else {
01251         jresp++;
01252       }
01253     }
01254     iresp++;
01255   }
01256 
01257   nvrt=GetVertexEntries();
01258   SortVertexByVsig();
01259   SortVertexByVsigInAllTrack();
01260 
01261   if (dblv > 0)
01262     cout <<nvrt<<" done"<<endl;
01263   PrintOutForDebug();
01264 
01265   //
01266   //now fit vertices and do chi**2 trimming
01267   //
01268  
01269   if (dblv > 0)
01270     cout << "   ---> Check 9: Chi**2 triming ... ";
01271   nvrt=GetVertexEntries();
01272   Int_t    ntracks_in_vertex =0;
01273   //Int_t    nvertices_in_track=0;
01274   Int_t    itrkmx            =0;
01275   Double_t maxchi2           =0.0;
01276   Double_t pchi2             =0.0;
01277   Int_t    itvrt             =0;
01278 
01279   ivrt = 0;
01280   Bool_t f_inc_iptrk;
01281   while (ivrt < nvrt) {
01282 
01283     vrt=GetVertexAt(ivrt);
01284     xvrt=vrt->GetVertexVector();
01285     vrt->VertexFit(xvrt);
01286 
01287     while(1) {
01288       ntracks_in_vertex=vrt->GetTrackEntries();
01289 
01290       f_inc_iptrk = (vrt->IndexOfTrack(iptrk) >= 0) ? kTRUE : kFALSE;
01291 
01292       if (ntracks_in_vertex < 2) {
01293 
01294         for (itrk=0 ; itrk < ntracks_in_vertex ; itrk++) {
01295           trk=vrt->GetTrackAt(itrk);
01296           // trk->RemoveVertex((TObject*)vrt);
01297           trk->RemoveVertex(vrt);
01298         }
01299         RemoveVertex(vrt);
01300         nvrt=GetVertexEntries();
01301         break;
01302 
01303       } else {
01304 
01305         //    xvinit = ipv;
01306         itrkmx =0;
01307         maxchi2=0.0;
01308         if (ntracks_in_vertex ==2 && f_inc_iptrk) {
01309           trk=0;
01310           for (itrk=0; itrk < ntracks_in_vertex ; itrk++) {
01311             trk=vrt->GetTrackAt(itrk);
01312             if (trk != iptrk) {
01313               break;
01314             }
01315           }
01316           xvrt=GetInitialVertexIP2T(vtrack_list->IndexOf(trk));
01317           vmaxi=GetVertexSignificance(xvrt);
01318           vrt->SetVertex(xvrt);
01319           vrt->SetVsig(vmaxi);
01320         } else {
01321           xvrt=vrt->GetVertexVector();
01322         }
01323         for (itrk=0; itrk < ntracks_in_vertex ; itrk++) {
01324           trk=vrt->GetTrackAt(itrk);
01325           if (trk != iptrk) {
01326             if ( (ntracks_in_vertex >= 3 && !f_inc_iptrk) ||
01327                  (ntracks_in_vertex >  3 &&  f_inc_iptrk)) {
01328               pchi2=vrt->GetPchi2(itrk);
01329             } else {
01330               pchi2=-2.0*
01331                 TMath::Log(TMath::Max(1e-40,
01332                                       trk->GetTrackProbabirity(xvrt,ipv)));
01333             }
01334           } else {
01335             if ( (ntracks_in_vertex >= 3 && !f_inc_iptrk) ||
01336                  (ntracks_in_vertex >  3 &&  f_inc_iptrk)) {
01337               pchi2=TMath::Power(vrt->GetDecayLength(ipv),2)
01338               /vrt->GetErrorDecayLength(ipv,iper);
01339             } else {
01340               pchi2=-2.0*TMath::Log(TMath::Max(1e-40,GetIPProbabirity(xvrt)));
01341             }
01342           }
01343           if (pchi2 > maxchi2) {
01344             maxchi2=pchi2;
01345             itrkmx =itrk;
01346           }
01347         }
01348 
01349         if (maxchi2 > xcut) {
01350           //Remove this vertex from track info.
01351           trk=vrt->GetTrackAt(itrkmx);
01352           // trk->RemoveVertex((TObject*)vrt);
01353           trk->RemoveVertex(vrt);
01354           //Drop this track
01355           vrt->RemoveTrackAt(itrkmx);
01356         } else {
01357           ivrt++;
01358           break;
01359         } // ifmaxchi2....
01360 
01361       } //more than 1 track in vertex
01362     } // loop over for chi**2 triming
01363   } //loop over vertices
01364 
01365   nvrt=GetVertexEntries();
01366   for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
01367     vrt=GetVertexAt(ivrt);
01368     vmaxi=GetVertexSignificance(vrt->GetVertexVector());
01369     vrt->SetVsig(vmaxi);
01370   }
01371   SortVertexByVsig();
01372   SortVertexByVsigInAllTrack();
01373 
01374   if (dblv > 0)
01375     cout <<nvrt<<" done"<<endl;
01376   PrintOutForDebug();
01377 
01378   //
01379   // now decide on track ambiguities by starting with highest spat. max in vrt
01380   //
01381   if (dblv > 0)
01382     cout << "   ---> Check 10 : Solve track ambiguities ... ";
01383   for (itrk=0 ; itrk < ntrk+1 ; itrk++ ) {
01384     //loop over tracks + ip
01385     trk  = GetTrackAt(itrk);
01386     vrti = 0;
01387     nvrt = trk->GetVertexEntries();
01388     for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
01389       vrt=(LCDVToplVRT*)(trk->GetVertexAt(ivrt));
01390       ntracks_in_vertex=vrt->GetTrackEntries();
01391       f_inc_iptrk = (vrt->IndexOfTrack(iptrk) >= 0) ? kTRUE : kFALSE;
01392       //require at least 2 trk in the vertex or in ip...
01393       if ( (ntracks_in_vertex > 1 && !f_inc_iptrk) || f_inc_iptrk) {
01394         vrti=vrt;
01395         break;
01396       }
01397     }
01398 
01399     itvrt=nvrt;
01400     while(--itvrt > 0) {
01401       vrtj=(LCDVToplVRT*)(trk->GetVertexAt(itvrt));
01402       if (vrtj == vrti) {
01403         vrtj = (LCDVToplVRT*)(trk->GetVertexAt(0));
01404         trk->AddVertexAt(vrtj,itvrt);
01405         trk->AddVertexAt(vrti,0    );
01406       }
01407       vrtj->RemoveTrack(trk,0);
01408       trk->RemoveVertexAt(itvrt);
01409     }
01410 
01411   }
01412 
01413   nvrt=GetVertexEntries();
01414   if (dblv > 0)
01415     cout<<nvrt<<" done"<<endl;
01416   PrintOutForDebug();
01417 
01418   //
01419   //re-sort vertices if at least two tracks or ip
01420   //
01421   if (dblv > 0)
01422     cout << "   ---> Check 11: Resort Vertices ...";
01423   nresp = 0;
01424   iresp = 0;
01425   ivrt=0;
01426   while(ivrt < nvrt) {
01427     vrt=GetVertexAt(ivrt);
01428     nvtrk=vrt->GetTrackEntries();
01429 
01430     f_inc_iptrk = (vrt->IndexOfTrack(iptrk) >= 0) ? kTRUE : kFALSE;
01431 
01432     //require at least 2 trk in the vertex or in ip...
01433     if ( (nvtrk > 1 && !f_inc_iptrk) || f_inc_iptrk) {
01434       vrt->VertexFit(ipv);
01435       ivrt++;
01436     } else {
01437       for (itrk=0 ; itrk < nvtrk ; itrk++) {
01438         trk=vrt->GetTrackAt(itrk);
01439         // trk->RemoveVertex((TObject*)vrt);
01440         trk->RemoveVertex(vrt);
01441       }
01442       RemoveVertex(vrt);
01443       nvrt=GetVertexEntries();
01444     }
01445   } // endloop of ivrt
01446 
01447   nvrt=GetVertexEntries();
01448   if (dblv > 0)
01449     cout <<nvrt<<" done"<<endl;
01450   PrintOutForDebug();
01451 
01452   // now find resolved vertices within rcut
01453   if (dblv > 0)
01454     cout << "   ---> Check 12: Resolve Vertices by Vsig (again)... ";
01455   nresp=GetVertexEntries();
01456   iresp=0;
01457   jresp=0;
01458   while(iresp < nresp-1) {
01459     vrti=GetVertexAt(iresp);
01460     jresp = iresp + 1;
01461     while(jresp < nresp) {
01462       vrtj=GetVertexAt(jresp);
01463       vrat=GetMinVrOverMinVs(vrti,vrtj,7);
01464 //      vrat=GetMinVrOverMinVs(vrti,vrtj,15);
01465       if (vrat > rcut) {
01466         MergeVertex(iresp,jresp,0);
01467         nresp=GetVertexEntries();
01468       } else {
01469         jresp++;
01470       }
01471     }
01472     iresp++;
01473   }
01474 
01475   Int_t nvrt_prev=nvrt;
01476   nvrt=GetVertexEntries();
01477   if (nvrt < nvrt_prev) {
01478     for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
01479       vrt=GetVertexAt(ivrt);
01480       xvrt=vrt->GetVertexVector();
01481       vrt->VertexFit(xvrt);    
01482     }
01483     SortVertexByVsig();
01484     //SortVertexByVsigInAllTrack();
01485   }
01486 
01487   if (dblv > 0)
01488     cout <<nvrt<<" done"<<endl;
01489   PrintOutForDebug();
01490 
01491   //
01492   //remove IP track...
01493   //
01494   if (dblv > 0)
01495     cout << "   ---> Check 13: Check IP ... ";
01496   Bool_t f_no_ipvrt=kTRUE;
01497   for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
01498     vrt=GetVertexAt(ivrt);
01499     ntracks_in_vertex=vrt->GetTrackEntries();
01500     //remove IP-track from vertices.
01501     itrk=0;
01502     while(itrk < ntracks_in_vertex) {
01503       trk=vrt->GetTrackAt(itrk);
01504       if (trk == iptrk) {
01505         f_no_ipvrt=kFALSE;
01506         vrt->RemoveTrackAt(itrk);
01507         ntracks_in_vertex=vrt->GetTrackEntries();
01508       } else {
01509         itrk++;
01510       }
01511     }
01512   }
01513   vtrack_list->Remove(iptrk);
01514   vtrack_list->Compress();
01515 
01516   //if (nvrt <= 0) {
01517   if (f_no_ipvrt || nvrt <= 0) {
01518     vrt=new((*vertex_list)[nvrt]) LCDVToplVRT(ipv);
01519     vmaxi=GetVertexSignificance(ipv);
01520     vrt->SetVsig(vmaxi);
01521     for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
01522       vrti=GetVertexAt(ivrt);
01523       vrat    =GetMinVrOverMinVs(vrt,vrti,7);
01524       if (vrat > rcut) {
01525         vrt=0;
01526         break;
01527       }
01528     }
01529     if (vrt == 0)
01530       RemoveVertex(vrt);
01531   }
01532 
01533   SortVertexByDecayLength();
01534   nvrt=GetVertexEntries();
01535   if (dblv > 0)
01536     cout <<nvrt<<" done"<<endl;
01537   PrintOutForDebug();
01538 
01539   //
01540   //save info.
01541   //
01542   if (dblv > 0)
01543     cout << "   ---> Check 14: Save info ... ";
01544   for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
01545     vrt=GetVertexAt(ivrt);
01546     vrt->CalcFourVector();
01547   }
01548 
01549   nvrt=GetVertexEntries();
01550   if (dblv > 0)
01551     cout <<nvrt<<" done"<<endl;
01552   if (nvrt == 0) { cout << " nvrt == 0 @ Check 12 " << ntrk << endl; }
01553 
01554   return 1;
01555 }
01556 //__________________________________________________________________________
01557 
01558 //-------------------------------------------
01559 // PT corrected mass relates...
01560 //-------------------------------------------
01561 
01562 //__________________________________________________________________________
01563 Int_t LCDVTopl::GetSeedVertex() {
01564   Int_t seedvertex_id=0;
01565   if (m_UseNN) {
01566     seedvertex_id=GetSeedVertexByNNMethod();
01567   } else {
01568     seedvertex_id=GetSeedVertexByClassicMethod();
01569   }
01570 
01571   return seedvertex_id;
01572 }
01573 //__________________________________________________________________________
01574 
01575 //__________________________________________________________________________
01576 Int_t LCDVTopl::GetSeedVertexByClassicMethod() {
01577   Int_t nseedvrt = -1;
01578   Int_t nvrt = GetVertexEntries();  
01579   if (nvrt < 2) return nseedvrt;
01580 
01581   Double_t maxvsig = -1.;
01582   LCDVToplVRT* vrt=GetVertexAt(0);
01583   for(Int_t ivrt =1; ivrt < nvrt ; ivrt++){
01584     vrt = GetVertexAt(ivrt);    
01585     
01586     // Apply cuts to select seed vertex candidates
01587 
01588     // Cut 1
01589     Double_t decayL  = vrt->GetDecayLength(ipv);
01590     if (decayL > m_SeedVrtRMax) continue;
01591 
01592     // Cut 2
01593     Int_t ntrks = vrt->GetTrackEntries();
01594     Double_t netCharge = vrt->GetVertexCharge();
01595     if (ntrks == 2 && netCharge == 0.) {      
01596       Double_t reconMass = vrt->GetMass();
01597       if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) continue;
01598     }
01599 
01600     // Cut 3
01601     TVector3 vrtpos(vrt->GetVertexVector() - ipv);
01602     if (vrtpos*pja < 0.) continue;
01603 
01604     // Cut 4
01605     TVector3 vrtmon(vrt->GetFourVector().Vect());
01606     if (vrtpos*vrtmon < 0.) continue;
01607 
01608     Double_t decayLerr = TMath::Sqrt(vrt->GetErrorDecayLength(ipv,iper));
01609     Double_t vrtsig = decayL/decayLerr;
01610     if (vrtsig < m_SeedVrtSigMin) continue;
01611 
01612     // Select the closest vertex
01613     Double_t vsig=vrt->GetVsig();
01614     if (nseedvrt < 0 || vsig > maxvsig) {
01615       maxvsig = vsig;  nseedvrt  = ivrt;
01616     }
01617   }  
01618 
01619   return nseedvrt;
01620 }
01621 //__________________________________________________________________________
01622 
01623 //__________________________________________________________________________
01624 Int_t LCDVTopl::GetSeedVertexByNNMethod() {
01625   Int_t nseedvrt = -1;
01626   Int_t nvrt = GetVertexEntries();  
01627   if (nvrt < 2) return nseedvrt;
01628 
01629   Double_t maxnn = -1.;
01630   LCDVToplVRT* vrt=GetVertexAt(0);
01631   Float_t nninp[3];
01632   Float_t nnout;
01633   for(Int_t ivrt =1; ivrt < nvrt ; ivrt++){
01634     vrt = GetVertexAt(ivrt);    
01635     
01636     // Apply cuts to select seed vertex candidates
01637 
01638     // Cut 1
01639     Double_t decayL  = vrt->GetDecayLength(ipv);
01640     TVector3 vrtpos(vrt->GetVertexVector() - ipv);
01641     TVector3 vrtmon(vrt->GetFourVector().Vect());
01642     Double_t decayLerr = TMath::Sqrt(vrt->GetErrorDecayLength(ipv,iper));
01643     Double_t vrtsig = decayL/decayLerr;
01644     Double_t pdang  = vrtpos.Angle(vrtmon);
01645     decayL=TMath::Log10(decayL);
01646     vrtsig=TMath::Log10(vrtsig);
01647     pdang =TMath::Log10(TMath::Max(pdang,1e-8));
01648     nninp[0]=TMath::Min(TMath::Max((decayL+3.0)/3.5,0.0),1.0);
01649     nninp[1]=TMath::Min(TMath::Max((vrtsig+0.5)/3.5,0.0),1.0);
01650     nninp[2]=TMath::Min(TMath::Max((pdang +3.5)/4.0,0.0),1.0);
01651     LCDNNSeedVertex(nninp,&nnout,0);
01652     if (nnout < m_SeedVrtNN) continue;
01653 
01654     // Cut 2
01655     Int_t ntrks = vrt->GetTrackEntries();
01656     Double_t netCharge = vrt->GetVertexCharge();
01657     if (ntrks == 2 && netCharge == 0.) {
01658       Double_t reconMass = vrt->GetMass();
01659       if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) continue;
01660     }
01661 
01662     // Select the closest vertex
01663     if (nseedvrt < 0 || nnout > maxnn) {
01664       maxnn = nnout;  nseedvrt  = ivrt;
01665     }
01666   }
01667 
01668   return nseedvrt;
01669 
01670 }
01671 //__________________________________________________________________________
01672 
01673 //__________________________________________________________________________
01674 void LCDVTopl::CalcPTMass(Int_t iseedvrt) {
01675   m_Rmass    = 0.;
01676   m_PtMass   = 0.;
01677   m_PtMassEr = 0.;
01678   m_Masstag  = 0.;
01679   m_MissPt   = 0.;
01680   
01681   if (iseedvrt < 0 ) { return; }
01682 
01683   TVector3 p3trk(0.0,0.0,0.0);
01684   const Double_t mpi = 0.139567;
01685 
01686   LCDVToplTRK* trk = 0;
01687   LCDVToplVRT* seedvrt = GetVertexAt(iseedvrt);
01688   Int_t nvtrk=seedvrt->GetTrackEntries();
01689   Int_t ivtrk=0;
01690   for (ivtrk=0 ; ivtrk < nvtrk ; ivtrk++) {
01691     trk = seedvrt->GetTrackAt(ivtrk);
01692     trk->SetTrkVal(1);
01693     AddVrtTrack(trk);
01694   }
01695   m_p4Vertex=seedvrt->GetFourVector();
01696 
01697   TVector3 xseedvrt(seedvrt->GetVertexVector());
01698 
01699   TVector3 ipseed(xseedvrt - ipv);
01700   Double_t dist = ipseed.Mag();
01701 
01702   // Attatch tracks to the Seed Vertex
01703   Float_t nninp[5];
01704   Float_t nnout=0;
01705   Int_t ntrks = GetTrackEntries();
01706   Double_t l   ;
01707   Double_t lodi;
01708   Double_t trdi;
01709   Double_t anta;
01710   Double_t ldist;
01711   for (Int_t itrk = 0; itrk < ntrks ; itrk++){
01712     trk = GetTrackAt(itrk);
01713 
01714     if (trk->GetTrkVal() != 0) continue;
01715 
01716     //Double_t l=trk->GetLminTrdiApprox(ipv,xseedvrt);
01717     //Double_t l=trk->GetLminTrdi(ipv,xseedvrt);
01718     //Double_t trdi = trk->GetTrdi(l,ipv,xseedvrt);
01719     //Double_t lodi = trk->GetLodi(l,ipv,xseedvrt);
01720 
01721     l    = 0;
01722     lodi = 0;
01723     trdi = 0;
01724     anta = 0;
01725     trk->GetTrdiLodiAnta(l,ipv,xseedvrt,&trdi,&lodi,&anta);
01726 
01727     ldist= lodi/dist;
01728 
01729     if (m_UseNN) {
01730       Double_t sbds = 
01731         TMath::Log10(trk->GetDistanceNormByError(l,ipv,iper));
01732       nninp[0]=
01733         TMath::Max(0.0,TMath::Min(1.0,trdi/0.3));
01734       nninp[1]=
01735         TMath::Max(0.0,TMath::Min(1.0,(lodi-0.3)/3.0));
01736       nninp[2]=
01737         TMath::Max(0.0,TMath::Min(1.0,(ldist-0.5)/3.0));
01738       nninp[3]=
01739         TMath::Max(0.0,TMath::Min(1.0,anta/TMath::Pi()));
01740       nninp[4]=
01741         TMath::Max(0.0,TMath::Min(1.0,sbds/3.0));
01742       LCDNNTrackAttach(nninp,&nnout,0);
01743       if (nnout < m_TrackAttNN) continue;
01744 
01745     } else {
01746       if (trdi  > m_TrdiMax ) continue;
01747       if (lodi  < m_LodiMin ) continue;
01748       if (ldist < m_LodrMin ) continue;
01749       if (ldist > m_LodrMax ) continue;
01750       if (anta  > m_AntaMax ) continue;
01751       if (trdi  > lodi      ) continue;
01752     }
01753     AddVrtTrack(trk);
01754     trk->SetTrkVal(2);
01755 
01756     p3trk = trk->GetMomentumVector(l); // Get momentum at Seed Vertex.
01757 
01758     m_p4Vertex += TLorentzVector(p3trk,TMath::Sqrt(p3trk.Mag2() + mpi*mpi));
01759 
01760     if (dblv > 9) {
01761       cout<<"      "<<" "
01762           <<"itrk="<<itrk<<" "
01763           <<"trdi="<<trdi<<" "
01764           <<"lodi="<<lodi<<" "
01765           <<"lovd="<<ldist<<" "
01766           <<endl;
01767     }
01768   }
01769 
01770   // Raw Mass
01771   m_Rmass = m_p4Vertex.Mag(); 
01772 
01773   // Momentum vector of reconstructed Vertex
01774   TVector3 p3Vertex=m_p4Vertex.Vect();
01775   Double_t pt = p3Vertex.Perp(ipseed);
01776 
01777   // Missing Pt
01778   m_MissPt = pt;
01779 
01780   // Pt corrected mass
01781   m_PtMass = TMath::Sqrt(m_Rmass*m_Rmass + pt*pt) + TMath::Abs(pt);
01782 
01783   //Find delta phi and delta lambda.
01784   Double_t sedxy = ipseed.Pt();
01785   Double_t pxy   = p3Vertex.Pt();
01786   Double_t dPhi  = p3Vertex.DeltaPhi(ipseed);
01787   Double_t dLam  = TMath::ATan((p3Vertex.Z()/pxy - ipseed.Z()/sedxy)/
01788                        (1.+p3Vertex.Z()*ipseed.Z()/(pxy*sedxy)));
01789 
01790   //Rotate vertex error matrix 
01791   TMatrixD evrt(3,3); seedvrt->GetErrorMatrixVertex(evrt);
01792   //TVector3 pvn   = p3Vertex;
01793   TVector3 pvn   = ipseed;
01794   Double_t pvnphi= pvn.Phi();
01795   Double_t csphi = TMath::Cos(pvnphi);
01796   Double_t snphi = TMath::Sin(pvnphi);
01797   Double_t csthe = pvn.CosTheta();
01798   Double_t snthe = TMath::Sqrt(1.0 - csthe*csthe);
01799   TMatrixD dxpdxT(3,3);
01800   dxpdxT(0,0)=+snphi; dxpdxT(0,1)=-csphi*csthe; dxpdxT(0,2)= csphi*snthe;
01801   dxpdxT(1,0)=-csphi; dxpdxT(1,1)=-snphi*csthe; dxpdxT(1,2)= snphi*snthe;
01802   dxpdxT(2,0)=+0.0  ; dxpdxT(2,1)=       snthe; dxpdxT(2,2)=       csthe;
01803   //dxpdxT(0,0)=+snphi; dxpdxT(1,0)=-csphi*csthe; dxpdxT(2,0)= csphi*snthe;
01804   //dxpdxT(0,1)=-csphi; dxpdxT(1,1)=-snphi*csthe; dxpdxT(2,1)= snphi*snthe;
01805   //dxpdxT(0,2)=+0.0  ; dxpdxT(1,2)=       snthe; dxpdxT(2,2)=       csthe;
01806   TMatrixD tmp0(evrt,TMatrixD::kMult,dxpdxT);
01807   TMatrixD verrot(tmp0,TMatrixD::kTransposeMult,dxpdxT);
01808 
01809   //Take errors into accound
01810   //Allow sigma IP error to make conservative Min missing P
01811   //Also allow seed vertex error n sigma Min missing pt
01812 
01813   Double_t dPhiM = TMath::Abs(dPhi) - 
01814     TMath::Sqrt(m_PTIP[0]*m_PTIP[0] + m_NVsigma*m_NVsigma*verrot(0,0)
01815          + m_NVsigma*m_NVsigma*verrot(2,2)
01816                 *TMath::Sin(dPhi)*TMath::Sin(dPhi) )/sedxy;
01817   if (dPhiM < 0.) dPhiM = 0.;
01818 
01819   Double_t dLamM = TMath::Abs(dLam) - 
01820     TMath::Sqrt(m_PTIP[1]*m_PTIP[1]*sedxy*sedxy/(dist*dist) 
01821                 + m_NVsigma*m_NVsigma*verrot(1,1)
01822                 + m_NVsigma*m_NVsigma
01823                 *verrot(2,2)*TMath::Sin(dLam)*TMath::Sin(dLam))/dist;
01824   if (dLamM < 0.) dLamM = 0.;
01825 
01826   pt = TMath::Sqrt(TMath::Power(p3Vertex.Mag()*TMath::Sin(dLamM),2)
01827                    +TMath::Power(pxy*TMath::Sin(dPhiM),2));
01828   
01829   // Ptcorrected mass with IP errors
01830   m_PtMassEr = TMath::Sqrt(m_Rmass*m_Rmass + pt*pt) + TMath::Abs(pt); 
01831 
01832   // Pt corrected IP err with 2*raw mass protection
01833   if (m_PtMassEr > 2.*m_Rmass) {
01834     m_Masstag = 2.*m_Rmass; 
01835   } else {
01836     m_Masstag = m_PtMassEr; 
01837   }
01838 }
01839 //__________________________________________________________________________
01840 
01841 //__________________________________________________________________________
01842 void LCDVTopl::CalcPTMass() {
01843   // Get Seed Vertex
01844   Int_t iseedvrt = GetSeedVertex();
01845   CalcPTMass(iseedvrt);
01846 }
01847 //__________________________________________________________________________
01848 
01849 //__________________________________________________________________________
01850 Double_t LCDVTopl::GetHFTagNN() {
01851   // Get Seed Vertex
01852   Int_t iseedvrt = GetSeedVertex();
01853   return GetHFTagNN(iseedvrt);
01854 }
01855 //__________________________________________________________________________
01856 
01857 //__________________________________________________________________________
01858 Double_t LCDVTopl::GetHFTagNN(Int_t iseed) {
01859   if (iseed < 0 ) { return -1.0; } //no seed vertex...
01860   if (m_Masstag == 0.0) {
01861     CalcPTMass(iseed);
01862   }
01863   LCDVToplVRT* seedvrt = GetVertexAt(iseed);
01864   TVector3     ipseed(seedvrt->GetVertexVector()-ipv);
01865   Double_t dist=ipseed.Mag();
01866   Double_t mvrt=GetMassTag();
01867   Double_t pvrt=15*mvrt-GetP4Vertex().Rho();
01868   Double_t dvrt=TMath::Log10(dist);
01869   Double_t tvrt=GetVrtTrackEntries();
01870   Double_t nvrt=GetVertexEntriesExceptV0();
01871   Float_t nninp[5];
01872   nninp[0]=TMath::Max(0.0,TMath::Min(1.0,mvrt/ 5.0));
01873   nninp[1]=TMath::Max(0.0,TMath::Min(1.0,(pvrt+20.0)/100.0));
01874   nninp[2]=TMath::Max(0.0,TMath::Min(1.0,(dvrt+2.0)/2.5));
01875   nninp[3]=TMath::Max(0.0,TMath::Min(1.0,tvrt/10.0));
01876   nninp[4]=TMath::Max(0.0,TMath::Min(1.0,(nvrt-1.0)/3.0));
01877   Float_t nnout;
01878   LCDNNHFSelect(nninp,&nnout,0);
01879   return nnout;
01880 }
01881 //__________________________________________________________________________
01882 
01883 //__________________________________________________________________________
01884 Int_t LCDVTopl::GetNsig(Double_t nsigma, TVector3 axis, TVector3 xfrom, 
01885                         TMatrixD& exfrom) {
01886   LCDVToplTRK* trk=0;
01887   Int_t itrk=0;
01888   Int_t ntrk=GetTrackEntries();
01889   Int_t nsig=0;
01890   Double_t l=0;
01891   for (itrk=0 ; itrk < ntrk ; itrk++) {
01892     trk=GetTrackAt(itrk);
01893     if (trk->GetV0flag()) continue;
01894     l=trk->CalcPOCA3(xfrom);
01895     if (trk->GetSignedDistanceNormByError(l,xfrom,axis,exfrom) > nsigma) {
01896       nsig++;
01897     }
01898   }
01899   return nsig;
01900 }
01901 //__________________________________________________________________________
01902 
01903 //__________________________________________________________________________
01904 Int_t LCDVTopl::GetNsig(Double_t nsigma) {
01905     return GetNsig(nsigma,pja,ipv,iper);
01906 }
01907 //__________________________________________________________________________
01908 
01909 //__________________________________________________________________________
01910 Int_t LCDVTopl::GetNsig(Double_t nsigma, TVector3 axis, TVector3 xfrom, 
01911                         TMatrixD& exfrom,TObjArray* nsig_list) {
01912   LCDVToplTRK* trk=0;
01913   Int_t itrk=0;
01914   Int_t ntrk=GetTrackEntries();
01915   Int_t nsig=0;
01916   Double_t l=0;
01917   for (itrk=0 ; itrk < ntrk ; itrk++) {
01918     trk=GetTrackAt(itrk);
01919     if (trk->GetV0flag()) continue;
01920     l=trk->CalcPOCA3(xfrom);
01921     if (trk->GetSignedDistanceNormByError(l,xfrom,axis,exfrom) > nsigma) {
01922       nsig_list->Add((TObject*)trk);
01923       nsig++;
01924     }
01925   }
01926   return nsig;
01927 }
01928 //__________________________________________________________________________
01929 
01930 //__________________________________________________________________________
01931 Int_t LCDVTopl::GetNsig(Double_t nsigma,TObjArray* nsig_list) {
01932     return GetNsig(nsigma,pja,ipv,iper,nsig_list);
01933 }
01934 //__________________________________________________________________________
01935 
01936 //__________________________________________________________________________
01937 Int_t LCDVTopl::GetNsig(Double_t nsigma, TVector3 axis, TVector3 xfrom, 
01938                         TMatrixD& exfrom,TObjArray* nsig_list,
01939                         Double_t rnsigcut, Double_t znsigcut) {
01940   LCDVToplTRK* trk=0;
01941   Int_t itrk=0;
01942   Int_t ntrk=GetTrackEntries();
01943   Int_t nsig=0;
01944   Double_t l=0;
01945   TVector3 ip2ca;
01946   Double_t rnsig;
01947   Double_t znsig;
01948   for (itrk=0 ; itrk < ntrk ; itrk++) {
01949     trk=GetTrackAt(itrk);
01950     if (trk->GetV0flag()) continue;
01951     //r & z check...
01952     ip2ca=trk->GetPositionVector(0.0)-xfrom;
01953     rnsig=ip2ca.Perp();
01954     znsig=TMath::Abs(ip2ca.Z());
01955     if (rnsig > rnsigcut) continue;
01956     if (znsig > znsigcut) continue;
01957     l=trk->CalcPOCA3(xfrom);
01958     if (trk->GetSignedDistanceNormByError(l,xfrom,axis,exfrom) > nsigma) {
01959       nsig_list->Add((TObject*)trk);
01960       nsig++;
01961     }
01962   }
01963   return nsig;
01964 }
01965 //__________________________________________________________________________
01966 
01967 //__________________________________________________________________________
01968 Int_t LCDVTopl::GetNsig(Double_t nsigma,TObjArray* nsig_list, 
01969                         Double_t rnsigcut, Double_t znsigcut) {
01970     return GetNsig(nsigma,pja,ipv,iper,nsig_list,rnsigcut,znsigcut);
01971 }
01972 //__________________________________________________________________________
01973 
01974 //__________________________________________________________________________
01975 void LCDVTopl::CheckV0Vertex() {
01976   Int_t ntrk=GetTrackEntries();
01977 
01978   LCDVToplTRK* trki=0;
01979   LCDVToplTRK* trkj=0;
01980   Int_t        itrk=0;
01981   Int_t        jtrk=0;
01982 
01983   TVector3 xvinit;
01984   LCDVToplVRT vrt;
01985   for (itrk=1 ; itrk < ntrk ; itrk++) { 
01986     trki=GetTrackAt(itrk);
01987     
01988     for (jtrk=0 ; jtrk < itrk ;jtrk++) { 
01989       trkj=GetTrackAt(jtrk);
01990 
01991       //track-track combination
01992       xvinit=GetInitialVertexT2T(itrk,jtrk);
01993 
01994       vrt.Clean();
01995       vrt.SetVertex(xvinit);
01996       vrt.AddTrack(trki);
01997       vrt.AddTrack(trkj);
01998 
01999       //V0 check...
02000       if (vrt.GetChi2() < 16.0) {
02001         vrt.CalcFourVector();
02002         if (CheckKs0Vertex(&vrt)) {
02003           trki->SetV0flag(1);
02004           trkj->SetV0flag(1);
02005         }
02006         if (CheckLambda0Vertex(&vrt)) {
02007           trki->SetV0flag(2);
02008           trkj->SetV0flag(2);
02009         }
02010       }
02011 
02012     }
02013   }
02014 }
02015 //__________________________________________________________________________
02016 
02017 //__________________________________________________________________________
02018 Bool_t LCDVTopl::CheckKs0Vertex(LCDVToplVRT* vrt) {
02019   Bool_t fk0s=kFALSE;
02020 
02021   Int_t    ntrks     = vrt->GetTrackEntries();
02022   Double_t netCharge = vrt->GetVertexCharge();
02023   if (ntrks == 2 && netCharge == 0.) {
02024     Double_t reconMass = vrt->GetMass();
02025     if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) fk0s=kTRUE;
02026   }
02027   
02028   return fk0s;
02029 }
02030 //__________________________________________________________________________
02031 
02032 //__________________________________________________________________________
02033 Bool_t LCDVTopl::CheckLambda0Vertex(LCDVToplVRT* vrt) {
02034   Bool_t flambda=kFALSE;
02035 
02036   Double_t mpi = 0.13957;
02037   Double_t mp  = 0.93827;
02038 
02039   Int_t    ntrks     = vrt->GetTrackEntries();
02040   Double_t netCharge = vrt->GetVertexCharge();
02041   if (ntrks == 2 && netCharge == 0.) {
02042     TVector3 p3trk1=vrt->GetTrackMomentumVector(0);
02043     TVector3 p3trk2=vrt->GetTrackMomentumVector(1);
02044     TLorentzVector p4trk1a(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mpi*mpi));
02045     TLorentzVector p4trk2a(p3trk2,TMath::Sqrt(p3trk1.Mag2()+mp *mp ));
02046     TLorentzVector p4vrta=p4trk1a+p4trk2a;
02047     if ( TMath::Abs(p4vrta.M() - 1.1157) < m_SeedVrtLambda0Mass ) {
02048       flambda=kTRUE;
02049     } else {
02050       TLorentzVector p4trk1b(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mp *mp ));
02051       TLorentzVector p4trk2b(p3trk2,TMath::Sqrt(p3trk1.Mag2()+mpi*mpi));
02052       TLorentzVector p4vrtb=p4trk1b+p4trk2b;
02053       if ( TMath::Abs(p4vrtb.M() - 1.1157) < m_SeedVrtLambda0Mass ) 
02054         flambda=kTRUE;
02055     }
02056   }
02057   
02058   return flambda;
02059 }
02060 //__________________________________________________________________________
02061 
02062 //__________________________________________________________________________
02063 void LCDVTopl::PrintOutForDebug() {
02064 
02065   if (dblv > 19) {
02066     Int_t nvrt=GetVertexEntries();
02067     LCDVToplVRT* vrt;
02068     Int_t        nvtrk;
02069     for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
02070       vrt=GetVertexAt(ivrt);
02071       nvtrk=vrt->GetTrackEntries();
02072       cout<<"    "
02073           <<"ivrt="<<ivrt<<" "
02074           <<"vrt="<<vrt<<" "
02075           <<"ntrk="<<nvtrk<<" "
02076           <<"decayL="<<vrt->GetDecayLength(ipv)<<" "
02077           <<"chi2="<<vrt->GetChi2()<<" "
02078           <<"Vsig="<<vrt->GetVsig()<<" "
02079           <<endl;
02080       if (dblv > 24) {
02081         cout<<"      "
02082             <<"vpos="
02083             <<vrt->GetVertexX()<<" "
02084             <<vrt->GetVertexY()<<" "
02085             <<vrt->GetVertexZ()<<" "
02086             <<endl;
02087       }
02088       if (dblv > 29) {
02089         for (Int_t itrk=0 ; itrk < nvtrk ; itrk++)
02090           cout<<"      "
02091               <<"itrk="<<itrk<<" "
02092               <<"trk="<<vrt->GetTrackAt(itrk)<<" "
02093               <<endl;
02094       }
02095     }
02096 
02097     if (dblv > 34) {
02098       LCDVToplTRK* trk;
02099       for (Int_t itrk=0 ; itrk < GetTrackEntries() ; itrk++) {
02100         trk=(LCDVToplTRK*)GetTrackAt(itrk);
02101         cout<<"      "
02102             <<"itrk="<<itrk<<" "
02103             <<"trk="<<trk<<" "
02104             <<endl;
02105         for (Int_t itvrt=0 ; itvrt < trk->GetVertexEntries() ; itvrt++) {
02106           vrt=(LCDVToplVRT*)(trk->GetVertexAt(itvrt));
02107           cout<<"    "
02108               <<"    "<<" "
02109               <<"itvrt="<<itvrt<<" "
02110               <<"vrt="<<vrt<<" "
02111               <<endl;
02112         }
02113       }
02114     }
02115 
02116   }
02117 
02118 }
02119 //__________________________________________________________________________
02120 

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