NcVertex.cxx

Go to the documentation of this file.

 
#include "NcVertex.h"
#include "Riostream.h"
 
ClassImp(NcVertex); // Class implementation to enable ROOT I/O
 
NcVertex::NcVertex() : NcJet(),NcPosition()
{
 
 Init();
 Reset();
 SetNvmax();
 SetNjmax();
}

void NcVertex::Init()
{
 
 fNvmax=0;
 fVertices=0;
 fConnects=0;
 fVertexCopy=0;
 fNjmax=0;
 fJets=0;
 fJetTracks=0;
 fJetCopy=0;
 fLines=0;
}

NcVertex::NcVertex(Int_t n) : NcJet(n),NcPosition()
{
 
 if (n<=0)
 {
  cout << " *** This NcJet initialisation was invoked via the NcVertex ctor." << endl;
 }
 Init();
 Reset();
 SetNvmax();
 SetNjmax();
}

NcVertex::~NcVertex()
{
 
 if (fVertices)
 {
  delete fVertices;
  fVertices=0;
 }
 if (fConnects)
 {
  delete fConnects;
  fConnects=0;
 }
 if (fJets)
 {
  delete fJets;
  fJets=0;
 }
 if (fJetTracks)
 {
  delete fJetTracks;
  fJetTracks=0;
 }
 if (fLines)
 {
  delete fLines;
  fLines=0;
 }
}

void NcVertex::SetOwner(Bool_t own)
{
 
 Int_t mode=1;
 if (!own) mode=0;
 if (fVertices) fVertices->SetOwner(own);
 fVertexCopy=mode;
 if (fJets) fJets->SetOwner(own);
 fJetCopy=mode;
 
 NcJet::SetOwner(own);
}

NcVertex::NcVertex(const NcVertex& v) : NcJet(v.fNtinit),NcPosition(v)
{
 
 Init();
 fNvtx=0;
 fNjets=0;
 SetNvmax(v.fNvmax);
 SetNjmax(v.fNjmax);
 SetTrackCopy(v.GetTrackCopy());
 SetVertexCopy(v.GetVertexCopy());
 SetJetCopy(v.GetJetCopy());
 SetId(v.GetId());
 
 // Copy all tracks except the ones coming from jets
 NcTrack* tx=0;
 Int_t jetflag=0,connect=0;
 NcTrack* tx2=0;
 for (Int_t it=1; it<=v.fNtrk; it++)
 {
  tx=v.GetTrack(it);
  if (tx)
  {
   jetflag=v.IsJetTrack(tx);
   connect=v.IsConnectTrack(tx);
 
   if (!jetflag && !connect) AddTrack(tx);
 
   if (connect)
   {
    if (!fConnects)
    {
     fConnects=new TObjArray(fNvmax);
     if (!fTrackCopy) fConnects->SetOwner();
    }
    tx2=new NcTrack(*tx);
    fConnects->Add(tx2);
    AddTrack(tx2,0);
   } 
  }
 }
 
 // Copy all the (secondary) vertices without re-creating connecting tracks
 // The connecting tracks have already been copied above
 NcVertex* vx=0;
 for (Int_t iv=1; iv<=v.GetNvertices(); iv++)
 {
  vx=v.GetVertex(iv);
  if (vx) AddVertex(vx,0); 
 }
 
 // Copy all the jets including the jet tracks for these jets for which
 // this was also the case in the original vertex
 NcJet* jx=0;
 for (Int_t ij=1; ij<=v.GetNjets(); ij++)
 {
  jx=v.GetJet(ij);
  if (jx)
  {
   jetflag=0;
   if (jx->GetNtracks())
   {
    tx=jx->GetTrack(1);
    if (tx)
    {
     jetflag=v.IsJetTrack(tx);
    }
   }
   AddJet(jx,jetflag);
  } 
 }
}

void NcVertex::SetNvmax(Int_t n)
{
 
 if (n > 0)
 {
  fNvmax=n;
 }
 else
 {
  fNvmax=1;
 }
 if (fVertices)
 {
  delete fVertices;
  fVertices=0;
 }
}

void NcVertex::SetNjmax(Int_t n)
{
 
 if (n > 0)
 {
  fNjmax=n;
 }
 else
 {
  fNjmax=1;
 }
 if (fJets)
 {
  delete fJets;
  fJets=0;
 }
}

void NcVertex::Reset()
{
 
 NcJet::Reset();
 
 ResetPosition();
 
 fNvtx=0;
 if (fNvmax>0) SetNvmax(fNvmax);
 if (fConnects)
 {
  delete fConnects;
  fConnects=0;
 }
 
 fNjets=0;
 if (fNjmax>0) SetNjmax(fNjmax);
 if (fJetTracks)
 {
  delete fJetTracks;
  fJetTracks=0;
 }
 
 if (fLines)
 {
  delete fLines;
  fLines=0;
 }
}

void NcVertex::ResetVertices()
{
 
 NcTrack* t;
 if (fConnects)
 {
  for (Int_t i=0; i<=fConnects->GetLast(); i++)
  {
   t=(NcTrack*)fConnects->At(i);
   NcTrack* test=(NcTrack*)fTracks->Remove(t);
   if (test)
   {
    fNtrk--;
    (Nc4Vector&)(*this)-=(Nc4Vector&)(*t);
    fQ-=t->GetCharge();
    if (fTrackCopy) delete t;
   }
  }
  fTracks->Compress();
 }
 
 fNvtx=0;
 if (fNvmax>0) SetNvmax(fNvmax);
 if (fConnects)
 {
  delete fConnects;
  fConnects=0;
 }
}

void NcVertex::AddJet(NcJet& j,Int_t tracks)
{
 
 if (!fJets)
 {
  fJets=new TObjArray(fNjmax);
  if (fJetCopy) fJets->SetOwner();
 }
 if (fNjets == fNjmax) // Check if maximum jet number is reached
 {
  fNjmax++;
  fJets->Expand(fNjmax);
 }
 
 // Add the jet to the list 
 NcJet* jx=&j;
 if (fJetCopy) jx=(NcJet*)j.Clone();
 
 if (jx)
 {
  fNjets++;
  fJets->Add(jx); 
 }
 
 // Add the tracks of the jet to this vertex
 if (tracks)
 {
  if (!fJetTracks)
  {
   fJetTracks=new TObjArray();
  }
  Int_t copy=1-(jx->GetTrackCopy());
  NcTrack* tj;
  for (Int_t i=1; i<=jx->GetNtracks(); i++)
  {
   tj=jx->GetTrack(i);
   if (tj)
   {
    AddTrack(tj,copy);
    fJetTracks->Add(tj);
   }
  }
 }
}

void NcVertex::AddVertex(NcVertex& v,Int_t connect)
{
 
 if (!fVertices)
 {
  fVertices=new TObjArray(fNvmax);
  if (fVertexCopy) fVertices->SetOwner();
 }
 if (fNvtx == fNvmax) // Check if maximum vertex number is reached
 {
  fNvmax++;
  fVertices->Expand(fNvmax);
 }
 
 // Add the linked (secondary) vertex to the list 
 NcVertex* vx=&v;
 if (fVertexCopy) vx=(NcVertex*)v.Clone();
 
 if (vx)
 {
  fNvtx++;
  fVertices->Add(vx);
 } 
 
 // Create connecting track and update 4-momentum and charge for current vertex
 if (connect)
 {
  NcTrack* t=new NcTrack();
  t->SetBeginPoint(GetPosition());
  t->SetEndPoint(v.GetPosition());
  t->SetCharge(v.GetCharge());
  t->Set4Momentum((Nc4Vector&)v);
 
  AddTrack(t,0);
 
  if (!fConnects)
  {
   fConnects=new TObjArray(fNvmax);
   if (!fTrackCopy) fConnects->SetOwner();
  }
  fConnects->Add(t);
 }
}

void NcVertex::Data(TString f,TString u)
{
 
 const char* name=GetName();
 const char* title=GetTitle();
 cout << " *NcVertex::Data*";
 if (strlen(name))  cout << " Name : " << GetName();
 if (strlen(title)) cout << " Title : " << GetTitle();
 cout << endl;
 cout << " Id : " << fUserId << " Invmass : " << GetInvmass()
      << " Charge : " << GetCharge() << " Momentum : " << GetMomentum()
      << " Ntracks : " << GetNtracks() << endl;
 cout << " Nvertices : " << fNvtx << " Njets : " << fNjets
      << " Energy scale : " << fEscale << " GeV" << endl;
 cout << " ";
 if (Nc4Vector::fV.HasVector()) Nc4Vector::Data(f,u);
 cout << "  Position";
 NcPosition::Data(f,u); 
} 

void NcVertex::List(TString f,TString u,TObjArray* tracks)
{
 
 Data(f,u); // Information of the current vertex
 
 Int_t ntk=0;
 if (tracks)
 {
  ntk=tracks->GetEntries();
 }
 else
 {
  ntk=GetNtracks();
 }
 
 if (!ntk)
 {
  cout << " *NcVertex::List* No tracks are present." << endl;
  return;
 }
 
 // The tracks of this vertex
 TObject* obj;
 NcTrack* tx; 
 for (Int_t it=1; it<=ntk; it++)
 {
  if (tracks)
  {
   obj=tracks->At(it-1);
   if (!obj) continue;
   if (!obj->InheritsFrom("NcTrack")) continue;
   tx=(NcTrack*)obj;
  }
  else
  {
   tx=GetTrack(it);
   if (!tx) continue;
  }
   cout << "  ---Track no. " << it << endl;
   cout << " ";
   tx->Data(f,u); 
 }
 
 // The secondary vertices of this vertex
 NcVertex* v; 
 for (Int_t iv=1; iv<=GetNvertices(); iv++)
 {
  v=GetVertex(iv);
  if (v)
  {
   cout << "  ---Level 1 sec. vertex no. " << iv << endl;
   cout << " ";
   v->Data(f,u); 
  }
  else
  {
   cout << " *NcVertex::List* Error : No sec. vertex present." << endl; 
  }
 }
} 

void NcVertex::ListAll(TString f,TString u,TObjArray* tracks)
{
 
 Data(f,u); // Information of the current vertex
 
 Int_t ntk=0;
 if (tracks)
 {
  ntk=tracks->GetEntries();
 }
 else
 {
  ntk=GetNtracks();
 }
 
 if (!ntk)
 {
  cout << " *NcVertex::ListAll* No tracks are present." << endl;
  return;
 }
 
 // The tracks of this vertex
 TObject* obj;
 NcTrack* tx; 
 for (Int_t it=1; it<=ntk; it++)
 {
  if (tracks)
  {
   obj=tracks->At(it-1);
   if (!obj) continue;
   if (!obj->InheritsFrom("NcTrack")) continue;
   tx=(NcTrack*)obj;
  }
  else
  {
   tx=GetTrack(it);
   if (!tx) continue;
  }
   cout << "  ---Track no. " << it << endl;
   cout << " ";
   tx->ListAll(f,u); 
 }
 
 NcVertex* v=this;
 Dumps(v,1,f,u); // Information of all sec. vertices
}

void NcVertex::Dumps(NcVertex* v,Int_t n,TString f,TString u)
{
 
 NcVertex* vs; 
 for (Int_t iv=1; iv<=v->GetNvertices(); iv++)
 {
  vs=v->GetVertex(iv);
  if (vs)
  {
   cout << "  ---Level " << n << " sec. vertex no. " << iv << endl;
   cout << " ";
   vs->Data(f,u); 
 
   // The tracks of this vertex
   NcTrack* t; 
   for (Int_t it=1; it<=vs->GetNtracks(); it++)
   {
    t=vs->GetTrack(it);
    if (t)
    {
     cout << "  ---Track no. " << it << endl;
     cout << " ";
     t->ListAll(f,u); 
    }
    else
    {
     cout << " *NcVertex::Dumps* Error : No track present." << endl; 
    }
   }
 
   // Go for next sec. vertex level of this sec. vertex recursively
   Dumps(vs,n+1,f,u);
  }
  else
  {
   cout << " *NcVertex::Dumps* Error : No sec. vertex present." << endl; 
  }
 }
} 

Int_t NcVertex::GetNvertices() const
{
 
 return fNvtx;
}

NcVertex* NcVertex::GetVertex(Int_t i) const
{
 
 if (!fVertices)
 {
  cout << " *NcVertex*::GetVertex* No (secondary) vertices present." << endl;
  return 0;
 }
 else
 {
  if (i<=0 || i>fNvtx)
  {
   cout << " *NcVertex*::GetVertex* Invalid argument i : " << i
        << " Nvtx = " << fNvtx << endl;
   return 0;
  }
  else
  {
   return (NcVertex*)fVertices->At(i-1);
  }
 }
}

NcVertex* NcVertex::GetIdVertex(Int_t id) const
{
 
 NcVertex* vx=0;
 NcVertex* v=0;
 if (!fVertices)
 {
  cout << " *NcVertex*::GetIdVertex* No (secondary) vertices present." << endl;
  return 0;
 }
 else
 {
  for (Int_t i=0; i<fNvtx; i++)
  {
   vx=(NcVertex*)fVertices->At(i);
   if (id == vx->GetId()) v=vx;
  }
  return v;
 }
}

void NcVertex::SetVertexCopy(Int_t j)
{
 
 if (!fVertices)
 {
  if (j==0 || j==1)
  {
   fVertexCopy=j;
  }
  else
  {
   cout << "*NcVertex::SetVertexCopy* Invalid argument : " << j << endl;
  }
 }
 else
 {
  cout << "*NcVertex::SetVertexCopy* Storage already contained vertices."
       << "  ==> VertexCopy mode not changed." << endl; 
 }
}

Int_t NcVertex::GetVertexCopy() const
{
 
 return fVertexCopy;
}

Int_t NcVertex::GetNjets() const
{
 
 return fNjets;
}

NcJet* NcVertex::GetJet(Int_t i) const
{
 
 if (!fJets)
 {
  cout << " *NcVertex*::GetJet* No jets present." << endl;
  return 0;
 }
 else
 {
  if (i<=0 || i>fNjets)
  {
   cout << " *NcVertex*::GetJet* Invalid argument i : " << i
        << " Njets = " << fNjets << endl;
   return 0;
  }
  else
  {
   return (NcJet*)fJets->At(i-1);
  }
 }
}

NcJet* NcVertex::GetIdJet(Int_t id) const
{
 
 NcJet* jx=0;
 NcJet* j=0;
 if (!fJets)
 {
  cout << " *NcVertex*::GetIdJet* No jets present." << endl;
  return 0;
 }
 else
 {
  for (Int_t i=0; i<fNjets; i++)
  {
   jx=(NcJet*)fJets->At(i);
   if (id == jx->GetId()) j=jx;
  }
  return j;
 }
}

void NcVertex::SetJetCopy(Int_t j)
{
 
 if (!fJets)
 {
  if (j==0 || j==1)
  {
   fJetCopy=j;
  }
  else
  {
   cout << "*NcVertex::SetJetCopy* Invalid argument : " << j << endl;
  }
 }
 else
 {
  cout << "*NcVertex::SetJetCopy* Storage already contained jets."
       << "  ==> JetCopy mode not changed." << endl; 
 }
}

Int_t NcVertex::GetJetCopy() const
{
 
 return fJetCopy;
}

Int_t NcVertex::IsConnectTrack(NcTrack* t) const
{
 
 Int_t connect=0;
 if (fConnects)
 {
  if (fConnects->FindObject(t)) connect=1;
 }
 return connect;
}

Int_t NcVertex::IsJetTrack(NcTrack* t) const
{
 
 Int_t jetflag=0;
 if (fJetTracks)
 {
  if (fJetTracks->FindObject(t)) jetflag=1;
 }
 return jetflag;
}

void NcVertex::Draw(Int_t secs,Int_t cons,Int_t jets)
{
 
 Double_t vec[3]={0,0,0};
 NcTrack* tx=0;
 NcVertex* vx=0;
 NcPosition* r=0;
 Nc3Vector p;
 Float_t charge;
 
 NcPosition dummy;
 
 if (fLines) delete fLines;
 fLines=new TObjArray();
 fLines->SetOwner();
 
 Int_t ntk=GetNtracks();
 for (Int_t jtk=1; jtk<=ntk; jtk++)
 {
  tx=GetTrack(jtk);
 
  if (!tx) continue;
 
  charge=tx->GetCharge();
 
  TPolyLine3D* line=new TPolyLine3D();
  fLines->Add(line);
 
  if (IsConnectTrack(tx))
  {
   if (cons==1)
   {
    r=tx->GetBeginPoint();
    if (!r) r=&dummy;
    r->GetPosition(vec,"car");
    line->SetNextPoint(vec[0],vec[1],vec[2]);
    r=tx->GetEndPoint();
    if (!r) r=&dummy;
    r->GetPosition(vec,"car");
    line->SetNextPoint(vec[0],vec[1],vec[2]);
    line->SetLineWidth(1);
   }
  }
  else
  {
   r=tx->GetClosestPoint();
   if (!r) r=&dummy;
   r->GetPosition(vec,"car");
   line->SetNextPoint(vec[0],vec[1],vec[2]);
   p=tx->Get3Momentum();
   p=p+(*r);
   p.GetVector(vec,"car");
   line->SetNextPoint(vec[0],vec[1],vec[2]);
   line->SetLineWidth(3);
  }
 
  line->SetLineColor(kGreen);              // Neutral track
  if (charge>0) line->SetLineColor(kRed);  // Positive track
  if (charge<0) line->SetLineColor(kBlue); // Negative track
 
  // Mark tracks belonging to jets
  if (IsJetTrack(tx))
  {
   if (jets==1)
   {
    line->SetLineWidth(2);
    line->SetLineColor(kMagenta);
   }
  }
    
  line->Draw();
 }
 
 // Go for secondary vertices if selected
 if (secs==1)
 {
  Int_t nvtx=GetNvertices();
  for (Int_t jvtx=1; jvtx<=nvtx; jvtx++)
  {
   vx=GetVertex(jvtx);
   if (vx) vx->Draw(secs,cons,jets);
  }
 }
}

TObjArray* NcVertex::SortJets(Int_t mode,TObjArray* jets,TObjArray* ordered)
{
 
 if (ordered) ordered->Clear();
 
 TObjArray* ajets=jets;
 if (!ajets) ajets=fJets;
 
 Int_t njets=0;
 if (ajets) njets=ajets->GetEntries();
 
 if (!mode || abs(mode)>13 || !jets || !njets) return 0;
 
 if (ordered)
 {
  ordered->Expand(njets);
 }
 else
 {
  if (fSelected) delete fSelected;
  fSelected=new TObjArray(njets);
 }
 
 Double_t val1,val2; // Values of the observable to be tested upon
 TObjArray* arr=0;   // Generic sorted array pointer to be used in the logic below
 
 if (ordered)
 {
  arr=ordered;
 }
 else
 {
  arr=fSelected;
 }
 
 TObject* obj=0;
 NcJet* jx=0;
 Int_t nord=0;
 for (Int_t i=0; i<njets; i++) // Loop over all jets of the array
 {
  obj=ajets->At(i);
  if (!obj) continue;
  if (!obj->InheritsFrom("NcJet")) continue;
  jx=(NcJet*)obj;
 
  if (nord == 0) // store the first jet at the first ordered position
  {
   nord++;
   arr->AddAt(jx,nord-1);
   continue;
  }
 
  for (Int_t j=0; j<=nord; j++) // put jet in the right ordered position
  {
   if (j == nord) // jet has smallest (mode<0) or largest (mode>0) observable value seen so far
   {
    nord++;
    arr->AddAt(jx,j); // add jet at the end
    break; // go for next jet
   }
 
   val1=0;
   val2=0;
   
   switch (abs(mode))
   {
    case 1:
     val1=jx->GetNtracks();
     val2=((NcJet*)arr->At(j))->GetNtracks();
     break;
    case 2:
     val1=jx->GetEnergy(1);
     val2=((NcJet*)arr->At(j))->GetEnergy(1);
     break;
    case 3:
     val1=jx->GetMomentum(1);
     val2=((NcJet*)arr->At(j))->GetMomentum(1);
     break;
    case 4:
     val1=jx->GetInvmass(1);
     val2=((NcJet*)arr->At(j))->GetInvmass(1);
     break;
    case 5:
     val1=jx->GetPt(1);
     val2=((NcJet*)arr->At(j))->GetPt(1);
     break;
    case 6:
     val1=jx->GetPl(1);
     val2=((NcJet*)arr->At(j))->GetPl(1);
     break;
    case 7:
     val1=jx->GetEt(1);
     val2=((NcJet*)arr->At(j))->GetEt(1);
     break;
    case 8:
     val1=jx->GetEl(1);
     val2=((NcJet*)arr->At(j))->GetEl(1);
     break;
    case 9:
     val1=jx->GetMt(1);
     val2=((NcJet*)arr->At(j))->GetMt(1);
     break;
    case 10:
     val1=jx->GetRapidity();
     val2=((NcJet*)arr->At(j))->GetRapidity();
     break;
    case 11:
     val1=jx->GetPseudoRapidity();
     val2=((NcJet*)arr->At(j))->GetPseudoRapidity();
     break;
    case 12:
     val1=jx->GetNsignals();
     val2=((NcJet*)arr->At(j))->GetNsignals();
     break;
    case 13:
     val1=jx->GetCharge();
     val2=((NcJet*)arr->At(j))->GetCharge();
     break;
   }
 
   if (mode<0 && val1 <= val2) continue;
   if (mode>0 && val1 >= val2) continue;
 
   nord++;
   for (Int_t k=nord-1; k>j; k--) // create empty position
   {
    arr->AddAt(arr->At(k-1),k);
   }
   arr->AddAt(jx,j); // put jet at empty position
   break; // go for next jet
  }
 }
 
 if (ordered)
 {
  return 0;
 }
 else
 {
  return fSelected;
 }
}

TObject* NcVertex::Clone(const char* name) const
{
 
 NcVertex* vtx=new NcVertex(*this);
 if (name)
 {
  if (strlen(name)) vtx->SetName(name);
 }
 return vtx;
}