////////////////////////////////////////////////////////////////
//
//  The TRootDSMC class is an extension of the RootDS class
//  that additionally allows to copy MC data into the 
//  mini-DSTs. For embedded data also the matching between
//  MC tracks and reconstructed tracks is stored.
//
///////////////////////////////////////////////////////////////

#include <iostream>
using namespace std;

// T49 includes
#include "TRootDSMC.h"
#include "T49VertexRoot.h"
#include "T49ParticleRoot.h"
#include "T49EventMCRoot.h"
#include "T49ParticleMCRoot.h"
#include "T49VertexMCRoot.h"

// ROOT includes
#include "TROOT.h"
#include "TTree.h"
#include "TList.h"
#include "TClonesArray.h"

// DSPACK includes
#include "dsc_prot.h"
#include "dspack_rootds.h"
#include "na49_event_mc_str_root.h"
#include "na49_event_str_root.h"
#include "xi_find_root.h"
#include "gteval_match_t_root.h"
#include "gteval_match_user_t_root.h"
#include "point_macros.h"      

#define MTL  0
#define MTR  1

ClassImp(TRootDSMC)

/*
$Log: TRootDSMC.C,v $
Revision 1.29  2005/05/20 11:24:59  cblume
Get rid of deprecated warnings

Revision 1.28  2003/12/18 16:48:11  cblume
Make the dynamic dimension of NSTrackSave proportional to all->n_vertex for Si+Si data

Revision 1.27  2003/09/18 15:18:31  cblume
Fix bug for MTV0 Xis

Revision 1.26  2003/09/17 14:09:11  cblume
Copy also the Xi-daughters from MTV0 correctly

Revision 1.25  2003/09/15 16:45:45  cblume
Fix bug in FillMCTrack()

Revision 1.24  2002/06/14 08:43:33  cblume
Change all fEvent to fEventMC

Revision 1.23  2002/05/06 07:56:15  cblume
Increase dimension of NSTrackSave

Revision 1.22  2002/02/04 16:09:52  cblume
Include T49VetoRoot information

Revision 1.21  2001/12/12 17:30:37  cblume
Change dimension of dynamic array. Add first and last point to MC tracks

Revision 1.20  2001/11/26 12:12:08  cblume
Change size of dynamic pointer arrays

Revision 1.19  2001/11/21 16:53:31  cblume
Set splitlevel to 99

Revision 1.18  2001/11/12 12:14:34  cblume
Update for ROOT v3.01

 * Revision 1.17  2000/10/27  12:26:33  cblume
 * Set the start vertex also in FillEvent() for more than four daughters
 *
 * Revision 1.16  2000/10/27  11:57:06  cblume
 * Change in verbose level for copying the WFA. Set start vertex properly when more than four daughters
 *
 * Revision 1.15  2000/10/27  09:00:04  cblume
 * Introduced switch for RFIO
 *
 * Revision 1.14  2000/10/25  09:33:56  cblume
 * Changed buffersize to 8000. Copy also WFA dataBug fix in matching of secondaries
 *
 * Revision 1.13  2000/10/15  01:00:57  cblume
 * Store also the number of matched points
 *
 * Revision 1.12  2000/10/09  10:25:10  cblume
 * General update
 *
 * Revision 1.11  2000/08/28  08:28:34  cblume
 * Adapted to new FillTrack()
 *
 * Revision 1.10  2000/04/17  13:57:42  cblume
 * Added fFillMCOnly and FillMCEvent()
 *
 * Revision 1.9  2000/04/14  14:26:56  cblume
 * Bug fix in FillEvent(). Changed TTree buffersize to 16000
 *
 * Revision 1.8  2000/04/13  15:19:51  cblume
 * Introduced cut on GSI V0s. Rtrack info for secondary tracks only 
 * when there is no corresponding main vertex track. 
 * FillRing() and FillBeam()
 *
 * Revision 1.7  2000/02/28  13:02:52  cblume
 * Introduced FillMCTrack(). Copy now also TOF information for MC tracks
 *
 * Revision 1.6  2000/02/23  15:42:38  cblume
 * Ok ok, lets allow four MC-daughter tracks
 *
 * Revision 1.5  2000/02/23  10:06:16  cblume
 * Added the possibility to store a third MC-daughter
 *
 * Revision 1.4  2000/02/08  14:04:39  cblume
 * Increased the dimensions of NSTrackSave and NVtxMCSave
 *
 * Revision 1.3  1999/12/16  12:24:59  cblume
 * Changed splitlevel in TTree::Branch() from 2 to 1
 *
 * Revision 1.2  1999/12/13  17:30:01  cblume
 * Cleanup, change in comments
 *
 * Revision 1.1  1999/11/23  12:49:01  cblume
 * Add src
 *
*/

//______________________________________________________________________________
 TRootDSMC::TRootDSMC():TRootDS() 
{
  //
  // TRootDSMC standard constructor
  //

  fFillMCTracks = kFALSE;
  fFillMCOnly   = kFALSE;
  
}

//______________________________________________________________________________
 TRootDSMC::~TRootDSMC() 
{
  //
  //  TRootDSMC destructor 
  //

  CloseDS();
    
}

//______________________________________________________________________________
 Bool_t TRootDSMC::ReadEventT49() 
{
  //
  //  Read the next event from the DSPACK input file
  //  and copies it directly into the T49 data structure.
  //  This allows to use the T49ANA classes on DSPACK data.
  //  Returns FALSE if end of input file is reached.
  //

  if (fEventMC) {
    if (fVerbose) printf("Clear ClonesArraysn");
    fEventMC->GetPrimaryParticles()->Delete(); 
    fEventMC->GetSecondaryParticles()->Delete(); 
    fEventMC->GetVertices()->Delete();
    fEventMC->GetSecondaryVertices()->Delete();
    fEventMC->GetMCParticles()->Delete();
    fEventMC->GetMCVertices()->Delete();
  }

  if (ReadEventDS()) {

    if (fVerbose) printf("Read one eventn");
    if (!(fEventMC)) fEventMC = new T49EventMCRoot(kTRUE);    

    if (fVerbose) printf("Fill one eventn");
    if (fFillMCOnly)
      return FillMCEvent();
    else
      return FillEvent();

  }
  else {

    return kFALSE;

  }

}

//______________________________________________________________________________
 Bool_t TRootDSMC::DS2Tree(char *InFname, char *OutFname, Int_t MaxEvents) 
{
  //
  //  Copy the events from the input file 'InFname' (DSPACK)
  //  to the ROOT mini-DST 'OutFname'. Up to MaxEvents are
  //  copied. The ROOT T49EventMCRoot objects are written to the 
  //  mini-DST in a ROOT TTree object, which allows access to 
  //  the event data in an n-tuple like fashion   
  //

  Char_t RootFname[256];
  Int_t  NumEvents = 0;   

  if (fRootFile == NULL) {

    fEventMC = new T49EventMCRoot(kTRUE);

    // Open the ROOT output via RFIO if enabled
    if (fRFIO) {
      sprintf(RootFname,"rfio:%s",OutFname);
    }
    else {
      sprintf(RootFname,"%s",OutFname);
    }    
    printf("Open File %sn",RootFname);
    fRootFile = TFile::Open(RootFname,"RECREATE","A Tree DST");

    printf("SetCompressionLevel(2)n");
    fRootFile->SetCompressionLevel(2);
    fRootFile->ls();
    printf("Create Treen");
    fTree = CreateTree();
    fTree->SetAutoSave(1990000000);
    if(fVerbose > 1) fTree->Print();

  }

  // Copy DS DST to Root file 
  cerr << "Open Input File:  " << InFname << endl;
  OpenDS(InFname);

  cerr << "Loop over events " << InFname << endl;
  while (GetDSpackOk() && NumEvents < MaxEvents) {

    if (ReadEventDS() == kTRUE) {
      if (fVerbose) printf("Fill one eventn");
      Bool_t eventFilled = kFALSE;
      if (fFillMCOnly)
        eventFilled = FillMCEvent();
      else
        eventFilled = FillEvent();
      if (eventFilled) {
        if (fVerbose) printf("Fill treen");
        fTree->Fill();
        if (fVerbose) printf("Clear ClonesArraysn");
        fEventMC->GetPrimaryParticles()->Delete(); 
        fEventMC->GetSecondaryParticles()->Delete(); 
        fEventMC->GetVertices()->Delete();
        fEventMC->GetSecondaryVertices()->Delete();
        fEventMC->GetMCParticles()->Delete();
        fEventMC->GetMCVertices()->Delete();
      }
      NumEvents++;
    }

  } 

  CloseDS();

  return kTRUE;

}

//______________________________________________________________________________
 TTree *TRootDSMC::CreateTree() 
{
  //
  //  Create the TTree object for the DS2Tree member
  //  function.
  //

  TTree   *tree;

  printf("Create NA49 Monte Carlo ROOT treen");
  tree = new TTree("T49MC","A NA49 Monte Carlo ROOT tree");

  if (fVerbose) {
    fEventMC->Dump();
    printf("Create MC event branchn");
  }
  //tree->Branch("Event","T49EventMCRoot",&fEventMC,64000,2);
  //tree->Branch("Event","T49EventMCRoot",&fEventMC,64000,1);
  //tree->Branch("Event","T49EventMCRoot",&fEventMC,16000,1);
  //tree->Branch("Event","T49EventMCRoot",&fEventMC,8000,1);
  // For ROOT v3.01/06
  tree->Branch("Event","T49EventMCRoot",&fEventMC,8000,99);
  if (fVerbose) {
    printf("tree okn");
  }

  return tree;

}

//______________________________________________________________________________
 Bool_t TRootDSMC::FillEvent() 
{
  //
  //  Copy event, track and point data into ROOT NA49
  //  objects for the current event.
  //
  //  Updated for V0-tracks and vertices (CBL 10/02/99)
  // 
  //  Included the Monte Carlo data (CBL 02/06/99)
  //

  Int_t          i;
  Int_t          Dummy;
  Int_t          TotalTracks = 0;
  Int_t          num;
  Int_t          num2;

  Char_t         name[256];

  event_t       *event_p;
  rtrack_t      *rtrack_p;
  track_t       *track_p;
  track_t       *track_tmp_p;
  vertex_t      *vtx_p;
  vertex_t      *vtx_list;
  vertex_t      *vtx_stop_p;
  vertex_t      *main_vtx;
  beam_t        *beam_p;
  ring_t        *ring_p;
  veto_t        *veto_p;
  pedestal_t    *veto_ped_p;
  tpc_t         *all_p;
  g_track_ref_t *track_ref = 0;
  g_ptr_list_t  *ptr_list  = 0;
  monitor_tpc_t *moni_p;
  bpd_vertex_t  *bpd_vtx_p;
  wfa_t         *wfa_p;
  wfa_time_t    *wfa_time_p;
    
  // For the Monte Carlo data    
  event_mc_t    *event_mc_p;
  etrack_mc_t   *etrack_mc_p;
  track_mc_t    *track_mc_p;
  vertex_mc_t   *vtx_mc_p;
  match_track_t *match_track_p;
  etrack_t      *etrack_p;

  static Int_t   nSkip = 0;

  T49ParticleRoot       *track;
  TClonesArray          *PrimaryTrackArray;
  TClonesArray          *SecondaryTrackArray;
  TObjArray             *VtxList;
  TObjArray             *SecVtxList;

  // For the Monte Carlo data
  T49ParticleMCRoot     *mcTrack;
  T49VertexMCRoot       *mcVertex;
  TClonesArray          *MCTrackArray;
  TObjArray             *MCVertexList;

  event_p = (event_t *) get_objects("event",&Dummy);
  if (event_p == NULL) {
    cerr << "No event_t in this event" << endl;
    return kFALSE;
  }

  event_mc_p = (event_mc_t *) get_objects("event_mc",&Dummy);
  if (event_mc_p == NULL) {
    cerr << "No event_mc_t in this event" << endl;
    return kFALSE;
  }

  all_p    = (tpc_t *) get_objects("all"   ,&num   );
  if (all_p == NULL) {
    cerr << "No all structure in this event" << endl;
    return kFALSE;
  }

  // Used to create the links between the different lists
  Int_t         NRTrackSave = TMath::Max( 250,all_p->n_rtrack + 1);
  Int_t         NSTrackSave = TMath::Max( 250,all_p->n_vertex * 3);
  Int_t         NVtxSave    = TMath::Max(5000,all_p->n_vertex * 3);
  Int_t         NVtxMCSave  = event_mc_p->n_vertex + 1;
  rtrack_t    **rtrack_save = new rtrack_t*[NRTrackSave];
  track_t     **strack_save = new track_t*[NSTrackSave];
  vertex_t    **vtx_save    = new vertex_t*[NVtxSave];
  vertex_mc_t **vtxMC_save  = new vertex_mc_t*[NVtxMCSave];
  for (i = 0; i < NRTrackSave; i++) rtrack_save[i] = NULL;
  for (i = 0; i < NSTrackSave; i++) strack_save[i] = NULL;
  for (i = 0; i < NVtxSave;    i++) vtx_save[i]    = NULL;
  for (i = 0; i < NVtxMCSave;  i++) vtxMC_save[i]  = NULL;
 
  if (nSkip > 0) {
    printf("Skip also the next eventn");
    nSkip--;
    return kFALSE;
  }

  // Skip problematic events
  for (Int_t iFound = 0; iFound < fSkipFound; iFound++) {
    if ((event_p->n_run   == fSkipRun[iFound]) &&
        (event_p->n_event == fSkipEvt[iFound])) {
      printf("Skip %d events starting with run %d, event %d.n"
            ,fSkipCnt[iFound],event_p->n_run,event_p->n_event);
      nSkip = fSkipCnt[iFound];
      nSkip--;
      return kFALSE;
    }
  }

  if (fVerbose) {
    printf("Fill event_t data (%d)n",fVerbose);
  }
  fEventMC->SetNRun(event_p->n_run);
  fEventMC->SetNEvent(event_p->n_event);
  fEventMC->SetDate(event_p->date);
  fEventMC->SetTime(event_p->time);
  fEventMC->SetTriggerMask(event_p->mask_trig);
  if (event_p->veto_p) {
    fEventMC->SetEveto(event_p->veto_p->eveto);
  }

  if (fVerbose) {
    printf("Copy monitor_tpc_t datan");
  }
  moni_p = (monitor_tpc_t *) get_objects("monitor_vt1",&Dummy);
  if (moni_p) {
    FillMonitor(fEventMC->GetMonitorTPC( 2),moni_p);
  }
  moni_p = (monitor_tpc_t *) get_objects("monitor_vt2",&Dummy);
  if (moni_p) {
    FillMonitor(fEventMC->GetMonitorTPC( 4),moni_p);
  }
  moni_p = (monitor_tpc_t *) get_objects("monitor_mtr",&Dummy);
  if (moni_p) {
    FillMonitor(fEventMC->GetMonitorTPC( 8),moni_p);
  }
  moni_p = (monitor_tpc_t *) get_objects("monitor_mtl",&Dummy);
  if (moni_p) {
    FillMonitor(fEventMC->GetMonitorTPC(16),moni_p);
  }

  if (fVerbose) {
    printf("Copy main vertex positionn");
  }
  bpd_vtx_p = (bpd_vertex_t *) get_objects("bpd_vertex",&Dummy);
  if (bpd_vtx_p) {
    if (fVerbose) {
      printf("Bpd_vertex structuren");
    }
    fEventMC->SetVertexX(bpd_vtx_p->x);
    fEventMC->SetVertexY(bpd_vtx_p->y);
    fEventMC->SetVertexZ(bpd_vtx_p->z);
    if (fVerbose) {
      printf("Create T49VertexRootn");
    }
    T49VertexRoot *Vertex = new T49VertexRoot(kFALSE);
    if (fVerbose) {
      printf("Fill vertex objectn");
    }
    Vertex->SetX(bpd_vtx_p->x);
    Vertex->SetY(bpd_vtx_p->y);
    Vertex->SetZ(bpd_vtx_p->z);
    Vertex->SetSigX(bpd_vtx_p->sx);
    Vertex->SetSigY(bpd_vtx_p->sy);
    Vertex->SetSigZ(bpd_vtx_p->sz);
    if (fVerbose) {
      printf("Add to vertex listn");
    }
    if ((VtxList = fEventMC->GetVertices())) {
      VtxList->Add(Vertex);
    }
  }
  else {
    printf("bpd_vertex = NULLn");
  }
  main_vtx = (vertex_t *) get_objects("vertex_fit",&Dummy);
  if (main_vtx) {
    if (fVerbose) {
      printf("Vertex_fit structuren");
    }
    fEventMC->SetVertexX(main_vtx->x);
    fEventMC->SetVertexY(main_vtx->y);
    fEventMC->SetVertexZ(main_vtx->z);
    if (fVerbose) {
      printf("Create T49VertexRootn");
    }
    T49VertexRoot *Vertex = new T49VertexRoot(kFALSE);
    if (fVerbose) {
      printf("Fill vertex objectn");
    }
    FillVertex(Vertex,main_vtx);
    if (fVerbose) {
      printf("Add to vertex listn");
    }
    if ((VtxList = fEventMC->GetVertices())) {
      VtxList->Add(Vertex);
    }
  }
  else {
    printf("vertex = NULLn");
  }

  if ((fEventMC->GetNEvent() % 100 == 0) || (fVerbose)) {
    sprintf(name,"event_%d_%d",fEventMC->GetNRun(),fEventMC->GetNEvent());
    sprintf(name,"Run %d, Event %d",fEventMC->GetNRun(),fEventMC->GetNEvent());
    printf("Event: %s - %g,%g,%gn",name
                                   ,fEventMC->GetVertexX()
                                   ,fEventMC->GetVertexY()
                                   ,fEventMC->GetVertexZ());
  }

  if (fVerbose) {
    printf("Copy the WFA datan");
  }
  wfa_p = (wfa_t *) get_objects("wfa",&Dummy);
  if (wfa_p) {
    fEventMC->SetWfaNbeam(wfa_p->n_beam);
    fEventMC->SetWfaNinter(wfa_p->n_inter);
  }
  else {
    if (fVerbose) {
      printf("No wfa_t foundn");
    }
  }
  wfa_time_p = (wfa_time_t *) get_objects("wfa_beam_time",&Dummy);
  if (wfa_time_p) {
    fEventMC->SetWfaBeamTime(wfa_time_p->wfa_time);
  }
  else {
    if (fVerbose) {
      printf("No wfa_beam_time foundn");
    }
  }
  wfa_time_p = (wfa_time_t *) get_objects("wfa_inter_time",&Dummy);
  if (wfa_time_p) {
    fEventMC->SetWfaInterTime(wfa_time_p->wfa_time);
  }
  else {
    if (fVerbose) {
      printf("No wfa_inter_time foundn");
    }
  }

  if (fVerbose) {
    printf("Copy beam_t datan");
  }
  beam_p = event_p->beam_p;
  if (!beam_p) beam_p = (beam_t *) get_objects("beam",&Dummy);
  if (beam_p && fEventMC->GetBeam()) {
    FillBeam(fEventMC->GetBeam(),beam_p);
  }
  else {
    printf("event_t->beam_p = %xn",(Int_t) event_p->beam_p);
  }

  if (fVerbose) {
    printf("Copy ring datan");
  }
  ring_p = event_p->ring_p;
  if (!ring_p) ring_p = (ring_t *) get_objects("ring",&Dummy);
  if (ring_p && fEventMC->GetRing()) {
    FillRing(fEventMC->GetRing(),ring_p);  
  }

  if (fVerbose) {
    printf("Copy veto datan");
  }
  veto_p = event_p->veto_p;
  if (!veto_p) veto_p = (veto_t *) get_objects("veto",&Dummy);
  veto_ped_p = (pedestal_t *) get_objects("veto_ped",&Dummy);
  if (veto_p && veto_ped_p && fEventMC->GetVeto()) {
    FillVeto(fEventMC->GetVeto(),veto_p,veto_ped_p);
  }

  if (fFillPrimary) {

    if (fVerbose) {
      printf("Copy primary track datan");
    }
    fNumTracks  = 0;
    TotalTracks = 0;

    track_ref = (g_track_ref_t *) get_objects("g_track_ref",&num);
    ptr_list  = (g_ptr_list_t  *) get_objects("g_ptr_list",&num2);

    if (!(PrimaryTrackArray = fEventMC->GetPrimaryParticles())) {
      printf("No primary TrackArray found in eventn");
    }

    if (fVerbose) {
      printf("Found g_track_ref_t at %x (%d tracks)n",(Int_t) track_ref,num);
    }
    if (num == num2) {

      for (i = 0; i < num; i++) {

        ptr_list = track_ref[i].g_ptr_list_p;
        if (!ptr_list) {
          fprintf(stderr,"ptr_list = NULL for %dn",i);
          break;
        }
        if (fVerbose > 1) {
          printf("%d: %xn",i,(Int_t) track_ref[i].g_ptr_list_p);
        }

        rtrack_p = ptr_list->rtrack_glb_p;
        if (fVerbose > 1) {
          printf("rtrack: %xn",(Int_t) rtrack_p);
        }
        if (!rtrack_p) {
          fprintf(stderr,"rtrack_p = NULL for %dn",i);
          break;
        }

        track_p = rtrack_p->first_track_p;
        if (fVerbose > 1) {
          printf("track: %xn",(Int_t) track_p);
        }
        if (!track_p) {
          fprintf(stderr,"track_p = NULL for %dn",i);
          break;
        }

        vtx_p = track_p->start_vertex_p;
        if (fVerbose > 1) {
          printf("vtx: %xn",(Int_t) vtx_p);
        }
        if (!vtx_p) {
          fprintf(stderr,"vtx_p = NULL for %dn",i);
          break;
        }           
        
        if ((track_p->pz           >    0) &&
            (rtrack_p->n_point     >   10) &&
            (track_p->iflag       ==    0) &&
            (track_p->rtrack_p    != NULL)) {

          // Copy tracks attached to the main vertex
          if (vtx_p->id_vtx == 0) {

            if (fVerbose > 1) {
              printf("Alloc T49ParticleRoot for main-vertex trackn");
            }
            track = new (PrimaryTrackArray->operator[](fNumTracks)) T49ParticleRoot();
            track->SetLabel(fNumTracks);
            if (fVerbose > 1) {
              printf("Call FillTrackn");
            }
            FillTrack(track,&track_ref[i],kTRUE);
            if (fVerbose > 1) {
              printf("FillTrack finishedn");
            }

            // Save the rtrack-pointer of this index for later
            if (fNumTracks < NRTrackSave) {
              rtrack_save[fNumTracks] = rtrack_p;
            }
            else {
              printf("Boundary error: fNumTracks = %d NRTrackSave = %dn"
                    ,fNumTracks,NRTrackSave);
  	    }

            fNumTracks++;

          }

        }

        TotalTracks++;

      }

    }
    else {
      printf("num = %d, num2 = %d, abort eventn",num,num2);
    }

    if (fVerbose) {
      printf("finished copying main-vertex tracksn");
    }

  }

  // Copy the secondary vertices and the corresponding daughter tracks
  Int_t  numVtx    = 0;
  Bool_t mainFound = kFALSE;
  if (fFillSecondary) {

    if (fVerbose) {
      if (fFillXis) {
        printf("Copy V0 and Xi datan");
      }
      else {
        printf("Copy V0 datan");
      }
    }
    // The TClonesArray for the non-vertex tracks
    if (!(SecondaryTrackArray = fEventMC->GetSecondaryParticles())) {
      printf("No secondary track array found in eventn");
    }
    // The TObjArray for the vertices
    if (!(SecVtxList = fEventMC->GetSecondaryVertices())) {
      printf("No secondary vertex list found in eventn");
    }

    fNumSecTracks   = 0;
    fNumSecVertices = 0;

    // New version of the loop to work around bug in xifind
    if ((all_p    = (tpc_t *)    get_objects("all"   ,&num   )) &&
        (vtx_list = (vertex_t *) get_objects("vertex",&numVtx)) &&
        (all_p->vertex_p)) {

      if (fFullSecondary) {
        track_ref = (g_track_ref_t *) get_objects("g_track_ref",&num);
        ptr_list  = (g_ptr_list_t  *) get_objects("g_ptr_list",&num2);
        if (num != num2) {
          printf("num = %d, num2 = %d, abort eventn",num,num2);
          return kFALSE;
	}
      }

      for (Int_t iVtx = all_p->vertex_p - vtx_list; iVtx < numVtx; iVtx++) {

        vtx_p = &vtx_list[iVtx];                                                  

	// Only good V0 vertices
        if (CheckV0(vtx_p)) {

          if (fVerbose > 1) {
            printf("Create T49VertexRootn");
          }
          T49VertexRoot *SecVertex = new T49VertexRoot(kFALSE);
          if (fVerbose > 1) {
            printf("Fill secondary vertex objectn");
          }
          FillSecVertex(SecVertex,vtx_p);
          if (fVerbose > 1) {
            printf("Add to vertex listn");
          }
          SecVtxList->Add(SecVertex);

          // Get the first daughter track
          track_p  = vtx_p->daughter_p;
          rtrack_p = track_p->rtrack_p;
          if (fVerbose > 1) {
            printf("Alloc T49ParticleRoot for first daughter trackn");
          }
          track = new (SecondaryTrackArray->operator[](fNumSecTracks)) T49ParticleRoot();
          track->SetLabel(fNumSecTracks);
	  // Save the index to the first particle in SecondaryTrackArray       
          SecVertex->SetSecondaryIndex(0,fNumSecTracks);
          if (fVerbose > 1) {
            printf("First particle index (non-vtx) = %dn",fNumSecTracks);
	  }

          // Save the index to the first particle in PrimaryTrackArray (main-vtx fit)       
          mainFound = kFALSE;
          for (i = 0; i < fNumTracks; i++) {
            if (rtrack_p == rtrack_save[i]) {
              SecVertex->SetPrimaryIndex(0,i);
              if (fVerbose > 1) {
                printf("First particle index (main-vtx) = %dn",i);
	      }
              mainFound = kTRUE;
              break;
	    }
          }

          if (fFullSecondary && !mainFound) {

            if (fVerbose > 1) {
              printf("Call FillTrack for first daughter trackn");
            }
            // Find the global track ref
            Bool_t track_ref_found = kFALSE;
            for (Int_t j = 0; j < num; j++) {
              ptr_list = track_ref[j].g_ptr_list_p;
              if (!ptr_list) {
                fprintf(stderr,"ptr_list = NULL for %dn",i);
                break;
              }
              if (rtrack_p == ptr_list->rtrack_glb_p) {
                if (!track_ref_found) {
                  FillTrack(track,&track_ref[j],kFALSE);
                  FillSecTrack(track,track_p,kTRUE);
                  track_ref_found = kTRUE;
                  if (fVerbose > 1) {
                    printf("FillTrack finishedn");
                  }
		}
                else {
                  fprintf(stderr,"More than one g_track_ref_t foundn");
		}
	      }
	    }
            if (!track_ref_found) {
              fprintf(stderr,"No g_track_ref_t foundn");
	    }

	  }
          else {

            if (fVerbose > 1) {
              printf("Call FillSecTrack for first daughter track. mainFound = %dn"
                    ,mainFound);
            }
            FillSecTrack(track, track_p, mainFound);
            if (fVerbose > 1) {
              printf("FillSecTrack finishedn");
            }

	  }

          if (rtrack_p->qpxz > 0) {
            if (fVerbose > 1) {
              printf("First track is positiven");
            }
            SecVertex->SetFirstTrackIsPositive(kTRUE);
	  }
	  else {
            if (fVerbose > 1) {
              printf("First track is negativen");
            }
            SecVertex->SetFirstTrackIsPositive(kFALSE);
          }
          // Save the track-pointer of this index for later
          if (fNumSecTracks < NSTrackSave) {
            strack_save[fNumSecTracks] = track_p;
	  }
          else {
            printf("Boundary error: fNumSecTracks = %d NSTrackSave = %dn"
                  ,fNumSecTracks,NSTrackSave);
	  }
          fNumSecTracks++;

          // Get the second daughter track
          track_p  = track_p->next_daughter_p;
          rtrack_p = track_p->rtrack_p;
          if (fVerbose > 1) {
            printf("Alloc T49ParticleRoot for second daughter trackn");
          }
          track = new (SecondaryTrackArray->operator[](fNumSecTracks)) T49ParticleRoot();
          track->SetLabel(fNumSecTracks);
	  // Save the index to the second particle in SecondaryTrackArray       
          SecVertex->SetSecondaryIndex(1,fNumSecTracks);
          if (fVerbose > 1) {
            printf("Second particle index (non-vtx) = %dn",fNumSecTracks);
	  }

          // Save the index to the second particle in PrimaryTrackArray (main-vtx fit)       
          mainFound = kFALSE;
          for (i = 0; i < fNumTracks; i++) {
            if (rtrack_p == rtrack_save[i]) {
              SecVertex->SetPrimaryIndex(1,i);
              if (fVerbose > 1) {
                printf("Second particle index (main-vtx) = %dn",i);
	      }
              mainFound = kTRUE;
              break;
	    }
          }

          if (fFullSecondary && !mainFound) {

            if (fVerbose > 1) {
              printf("Call FillTrack for second daughter trackn");
            }
            // Find the global track ref
            Bool_t track_ref_found = kFALSE;
            for (Int_t j = 0; j < num; j++) {
              ptr_list = track_ref[j].g_ptr_list_p;
              if (!ptr_list) {
                fprintf(stderr,"ptr_list = NULL for %dn",i);
                break;
              }
              if (rtrack_p == ptr_list->rtrack_glb_p) {
                if (!track_ref_found) {
                  FillTrack(track,&track_ref[j],kFALSE);
                  FillSecTrack(track,track_p,kTRUE);
                  track_ref_found = kTRUE;
                  if (fVerbose > 1) {
                    printf("FillTrack finishedn");
                  }
		}
                else {
                  fprintf(stderr,"More than one g_track_ref_t foundn");
		}
	      }
	    }
            if (!track_ref_found) {
              fprintf(stderr,"No g_track_ref_t foundn");
	    }

	  }
          else {

            if (fVerbose > 1) {
              printf("Call FillSecTrack for second daughter track. mainFound = %dn"
                    ,mainFound);
            }
            FillSecTrack(track,track_p,mainFound);
            if (fVerbose > 1) {
              printf("FillSecTrack finishedn");
            }

	  }
          // Save the track-pointer of this index for later
          if (fNumSecTracks < NSTrackSave) {
            strack_save[fNumSecTracks] = track_p;
	  }
          else {
            printf("Boundary error: fNumSecTracks = %d NSTrackSave = %dn"
                  ,fNumSecTracks,NSTrackSave);
	  }
          fNumSecTracks++;

          fNumSecVertices++;            

	}

        // The Xis
        if ((fFillXis) && (CheckXi(vtx_p))) {

          if (fVerbose > 1) {
            printf("Create T49VertexRootn");
          }
          T49VertexRoot *SecVertex = new T49VertexRoot(kFALSE);
          if (fVerbose > 1) {
            printf("Fill secondary vertex object with Xin");
          }
          FillXiVertex(SecVertex,vtx_p);
          if (fVerbose > 1) {
            printf("Add to vertex listn");
          }
          SecVtxList->Add(SecVertex);
          if (fNumSecVertices < NVtxSave) {
            vtx_save[fNumSecVertices] = vtx_p;
          }

          // Get the first daughter track
          track_p  = vtx_p->daughter_p;
	  // If Xi was found by MTV0 the charged track is the second one
          if (vtx_p->iflag & 0x2000) {
            track_p = track_p->next_daughter_p;
            if (!track_p) {
              printf("Error: No charged track from MTV0-Xin");
              continue;
	    }
	  }
          rtrack_p = track_p->rtrack_p;

          // Reconstruct the linked list for the tracks produced by the xifinder
          track_tmp_p = rtrack_p->first_track_p;
	  while (track_tmp_p->next_track_p) {
            track_tmp_p = track_tmp_p->next_track_p;
	  }
          track_tmp_p->next_track_p = track_p;

          if (fVerbose > 1) {
            printf("Alloc T49ParticleRoot for first Xi daughter trackn");
          }
          track = new (SecondaryTrackArray->operator[](fNumSecTracks)) T49ParticleRoot();
          track->SetLabel(fNumSecTracks);
	  // Save the index to the first particle in SecondaryTrackArray       
          SecVertex->SetSecondaryIndex(0,fNumSecTracks);
          if (fVerbose > 1) {
            printf("First particle index (non-vtx) = %dn",fNumSecTracks);
	  }

          // Save the index to the first particle in PrimaryTrackArray (main-vtx fit)       
          mainFound = kFALSE;
          if (rtrack_p) {
            for (i = 0; i < fNumTracks; i++) {
              if (rtrack_p == rtrack_save[i]) {
                SecVertex->SetPrimaryIndex(0,i);
                if (fVerbose > 1) {
                  printf("First particle index (main-vtx) = %dn",i);
	        }
                mainFound = kTRUE;
                break;
	      }
	    }
          }

          if (fFullSecondary && !mainFound && rtrack_p) {

            if (fVerbose > 1) {
              printf("Call FillTrack for first Xi daughter trackn");
            }
            // Find the global track ref
            Bool_t track_ref_found = kFALSE;
            for (Int_t j = 0; j < num; j++) {
              ptr_list = track_ref[j].g_ptr_list_p;
              if (!ptr_list) {
                fprintf(stderr,"ptr_list = NULL for %dn",i);
                break;
              }
              if (rtrack_p == ptr_list->rtrack_glb_p) {
                if (!track_ref_found) {
                  FillTrack(track,&track_ref[j],kFALSE);
                  FillSecTrack(track,track_p,kTRUE);
                  track_ref_found = kTRUE;
                  if (fVerbose > 1) {
                    printf("FillTrack finishedn");
                  }
		}
                else {
                  fprintf(stderr,"More than one g_track_ref_t foundn");
		}
	      }
	    }
            if (!track_ref_found) {
              fprintf(stderr,"No g_track_ref_t foundn");
	    }

	  }
          else {

            if (fVerbose > 1) {
              printf("Call FillSecTrack for first Xi daughter track. mainFound = %dn"
                    ,mainFound);
            }
            FillSecTrack(track, track_p, mainFound);
            if (fVerbose > 1) {
              printf("FillSecTrack finishedn");
            }

	  }

          if (track_p->qpxz > 0) {
            if (fVerbose > 1) {
              printf("First track is positiven");
            }
            SecVertex->SetFirstTrackIsPositive(kTRUE);
	  }
	  else {
            if (fVerbose > 1) {
              printf("First track is negativen");
            }
            SecVertex->SetFirstTrackIsPositive(kFALSE);
          }
          // Save the track-pointer of this index for later
          if (fNumSecTracks < NSTrackSave) {
            strack_save[fNumSecTracks] = track_p;
	  }
          else {
            printf("Boundary error: fNumSecTracks = %d NSTrackSave = %dn"
                  ,fNumSecTracks,NSTrackSave);
	  }
          fNumSecTracks++;

          // Get the second daughter track (neutral track / V0)
	  // If Xi was found by MTV0 the neutral track is the first one
          if (vtx_p->iflag & 0x2000) {
            track_p  = vtx_p->daughter_p;
	  }
          else {
            track_p  = track_p->next_daughter_p;
	  }
          rtrack_p = track_p->rtrack_p;
          if (fVerbose > 1) {
            printf("Alloc T49ParticleRoot for second Xi daughter trackn");
          }
          track = new (SecondaryTrackArray->operator[](fNumSecTracks)) T49ParticleRoot();
          track->SetLabel(fNumSecTracks);
	  // Save the index to the second particle in SecondaryTrackArray       
          SecVertex->SetSecondaryIndex(1,fNumSecTracks);
          if (fVerbose > 1) {
            printf("Second particle index (non-vtx) = %dn",fNumSecTracks);
	  }

	  // No rtrack information for the Lambda track 
          mainFound = kFALSE;
          if (fVerbose > 1) {
            printf("Call FillSecTrack for second Xi daughter track. mainFound = %dn"
                  ,mainFound);
          }
          FillSecTrack(track,track_p,mainFound);
          if (fVerbose > 1) {
            printf("FillSecTrack finishedn");
          }
          // Save the track-pointer of this index for later
          if (fNumSecTracks < NSTrackSave) {
            strack_save[fNumSecTracks] = track_p;
	  }
          else {
            printf("Boundary error: fNumSecTracks = %d NSTrackSave = %dn"
                  ,fNumSecTracks,NSTrackSave);
	  }
          fNumSecTracks++;

          fNumSecVertices++;            

	  // Find the Lambda vertex and add it to the list if neccessary
          vtx_stop_p = track_p->stop_vertex_p;
          if (vtx_stop_p) {
            // Check whether this V0 is already copied
            Bool_t foundV0 = kFALSE;
            Int_t  indexV0 = 0;
            for (i = 0; i < fNumSecVertices; i++) {
              if (vtx_stop_p == vtx_save[i]) {
                foundV0 = kTRUE;
                indexV0 = i;
                break;
	      }
            }
            if (foundV0) {
              track->SetStopVertexIndex(indexV0);
	    }
            else {
              track->SetStopVertexIndex(fNumSecVertices);
              // Add the V0 to the secondary vertex list
              T49VertexRoot *SecVertex = new T49VertexRoot(kFALSE);
              FillSecVertex(SecVertex,vtx_stop_p);
              track_p  = vtx_stop_p->daughter_p;
              rtrack_p = track_p->rtrack_p;
              track_tmp_p = rtrack_p->first_track_p;
	      while (track_tmp_p->next_track_p) {
                track_tmp_p = track_tmp_p->next_track_p;
	      }
              track_tmp_p->next_track_p = track_p;
              track = new (SecondaryTrackArray->operator[](fNumSecTracks)) T49ParticleRoot();
              track->SetLabel(fNumSecTracks);
              SecVertex->SetSecondaryIndex(0,fNumSecTracks);
              Bool_t track_ref_found = kFALSE;
              for (Int_t j = 0; j < num; j++) {
                ptr_list = track_ref[j].g_ptr_list_p;
                if (rtrack_p == ptr_list->rtrack_glb_p) {
                  if (!track_ref_found) {
                    FillTrack(track,&track_ref[j],kFALSE);
                    FillSecTrack(track,track_p,kTRUE);
                    track_ref_found = kTRUE;
		  }
	        }
	      }
              if (rtrack_p->qpxz > 0) {
                SecVertex->SetFirstTrackIsPositive(kTRUE);
	      }
	      else {
                SecVertex->SetFirstTrackIsPositive(kFALSE);
              }
              if (fNumSecTracks < NSTrackSave) {
                strack_save[fNumSecTracks] = track_p;
	      }
              else {
                printf("Boundary error: fNumSecTracks = %d NSTrackSave = %dn"
                      ,fNumSecTracks,NSTrackSave);
	      }
              fNumSecTracks++;
              track_p  = track_p->next_daughter_p;
              rtrack_p = track_p->rtrack_p;
              track_tmp_p = rtrack_p->first_track_p;
	      while (track_tmp_p->next_track_p) {
                track_tmp_p = track_tmp_p->next_track_p;
	      }
              track_tmp_p->next_track_p = track_p;
              track = new (SecondaryTrackArray->operator[](fNumSecTracks)) T49ParticleRoot();
              track->SetLabel(fNumSecTracks);
              SecVertex->SetSecondaryIndex(1,fNumSecTracks);
              track_ref_found = kFALSE;
              for (Int_t j = 0; j < num; j++) {
                ptr_list = track_ref[j].g_ptr_list_p;
                if (rtrack_p == ptr_list->rtrack_glb_p) {
                  if (!track_ref_found) {
                    FillTrack(track,&track_ref[j],kFALSE);
                    FillSecTrack(track,track_p,kTRUE);
                    track_ref_found = kTRUE;
		  }
	        }
	      }
              if (fNumSecTracks < NSTrackSave) {
                strack_save[fNumSecTracks] = track_p;
	      }
              else {
                printf("Boundary error: fNumSecTracks = %d NSTrackSave = %dn"
                      ,fNumSecTracks,NSTrackSave);
	      }
              fNumSecTracks++;
              SecVtxList->Add(SecVertex);
              if (fNumSecVertices < NVtxSave) {
                vtx_save[fNumSecVertices] = vtx_stop_p;
              }
              fNumSecVertices++;            
	    }

	  }
          else {
            fprintf(stderr,"No stop vertex foundn");
	  } 

	}

        //vtx_p = vtx_p->next_p;
        //vtx_p++;

      }

    }
    else {
      printf("No vertices foundn");
    }
      
    if (fVerbose) {
      printf("Finished copying V0 datan");
    }

  }

  // Copy MC data with the links to the matching reconstructed tracks
  if (fFillMCTracks) {

    fNumMCTracks   = 0;
    fNumMCVertices = 0;

    if (fVerbose) {
      printf("Copy Monte Carlo datan");
    }
    // The TClonesArray for the MC tracks
    if (!(MCTrackArray  = fEventMC->GetMCParticles())) {
      printf("No Monte Carlo track array found in eventn");
    }
    // The TObjArray for the MC vertices
    if (!(MCVertexList  = fEventMC->GetMCVertices())) {
      printf("No Monte Carlo vertex list found in eventn");
    }

    // Get pointer to the etrack_mc_t structure 
    etrack_mc_p = (etrack_mc_t *) get_objects("etrack_mc",&num);

    // Loop over all MC tracks.
    for (Int_t inum = 0; inum < num; inum++) {

      // Get pointer to the MC track
      track_mc_p = etrack_mc_p[inum].track_p;

      // Create a new MC track and add it to the track array
      if (fVerbose > 1) printf("Create T49ParticleMCRoot for pid %dn",track_mc_p->pid);
      mcTrack = new (MCTrackArray->operator[](fNumMCTracks)) T49ParticleMCRoot();	

      // Fill the MC track
      FillMCTrack(mcTrack,track_mc_p);

      if (fVerbose > 1) {
        printf("match_track_p = %xn",(Int_t) etrack_mc_p[inum].match_p);
      }

      // Set the links to the matching reconstructed tracks
      Int_t iMatchPri = 0;
      Int_t iMatchSec = 0;
      // Loop through all matched global rtracks
      Int_t iMatch    = 0;
      for (match_track_p = etrack_mc_p[inum].match_p;
           match_track_p;
           match_track_p = match_track_p->next_match_p) {
        
        iMatch++;
        if (iMatch > match_track_p->n_match) break;

        if ((etrack_p  = match_track_p->etrack_p  ) &&
            (track_ref = etrack_p->gtrack_p       ) &&
            (ptr_list  = track_ref->g_ptr_list_p  ) &&
            (rtrack_p  = ptr_list->rtrack_glb_p   )) {

	  // Loop through all tracks associated to this rtrack
          for (track_p = rtrack_p->first_track_p;
               track_p;
               track_p = track_p->next_track_p) {
            vtx_p = track_p->start_vertex_p;
            if (vtx_p->id_vtx == 0) {
              if (fVerbose > 1) printf("Looking for matched primary trackn");
              for (i = 0; i < fNumTracks; i++) {
                if (rtrack_p == rtrack_save[i]) {
                  //mcTrack->SetPriMatched(i,etrack_p->n_pnt_match);
                  mcTrack->SetPriMatched(i,match_track_p->n_pnt_common);
                  iMatchPri++;
                  if (fVerbose > 1)
                    printf("Found matched primary track. Index = %dn",i);
		}
	      }
	    }
            else {
              if (fVerbose > 1) printf("Looking for matched secondary trackn");
              for (i = 0; i < fNumSecTracks; i++) {
                if (track_p == strack_save[i]) {
                  //mcTrack->SetSecMatched(i,etrack_p->n_pnt_match);
                  mcTrack->SetSecMatched(i,match_track_p->n_pnt_common);
                  iMatchSec++;
                  if (fVerbose > 1)
                    printf("Found matched secondary track. Index = %dn",i);
	        }
	      }
	    }
	  }

	}

      }
      if (fVerbose > 1) {
        printf("Found %d matched primary tracksn"  ,iMatchPri);
        printf("Found %d matched secondary tracksn",iMatchSec);  
      }

      // Save the start-vertex
      if (fVerbose > 1) printf("Save the start vertexn");
      if ((vtx_mc_p = track_mc_p->start_vertex_p)) {
        Int_t indexVtx = -1;
        for (i = 0; i < fNumMCVertices; i++) {
          if (vtx_mc_p == vtxMC_save[i]) {
            indexVtx = i;
            break;
	  }
        }
        if (indexVtx >= 0) {
          mcVertex = (T49VertexMCRoot *) MCVertexList->At(indexVtx);
          if (fVerbose > 1) printf("Vertex already copied %dn",indexVtx);
        }
        else {
          if (fVerbose > 1) printf("Create new T49VertexMCRootn");
          mcVertex = new T49VertexMCRoot();
          MCVertexList->Add(mcVertex);        
          mcVertex->SetX(vtx_mc_p->x);
          mcVertex->SetY(vtx_mc_p->y);
          mcVertex->SetZ(vtx_mc_p->z);
          vtxMC_save[fNumMCVertices] = vtx_mc_p;
          indexVtx = fNumMCVertices;
          fNumMCVertices++;
        }
        if      (mcVertex->GetFirstDaughter()  < 0) {
          if (fVerbose > 1) printf("Set first daughtern");
          mcVertex->SetFirstDaughter(fNumMCTracks);
          mcTrack->SetStartVertex(indexVtx);
        }
        else if (mcVertex->GetSecondDaughter() < 0) {
          if (fVerbose > 1) printf("Set second daughtern");
          mcVertex->SetSecondDaughter(fNumMCTracks);
          mcTrack->SetStartVertex(indexVtx);
        }
        else if (mcVertex->GetThirdDaughter()  < 0) {
          if (fVerbose > 1) printf("Set third daughtern");
          mcVertex->SetThirdDaughter(fNumMCTracks);
          mcTrack->SetStartVertex(indexVtx);
        }
        else if (mcVertex->GetForthDaughter()  < 0) {
          if (fVerbose > 1) printf("Set forth daughtern");
          mcVertex->SetForthDaughter(fNumMCTracks);
          mcTrack->SetStartVertex(indexVtx);
        }
        else {
          if (fVerbose) printf("More than four MC-daughters for MC-vertex %dn"
                              ,indexVtx);
          mcTrack->SetStartVertex(indexVtx);
        }
      }
      else {
        if (fVerbose > 1) printf("No start vertex foundn");
      }

      // Save the stop vertex
      if (fVerbose > 1) printf("Save the stop vertexn");
      if ((vtx_mc_p = track_mc_p->stop_vertex_p)) {
        Int_t indexVtx = -1;
        for (i = 0; i < fNumMCVertices; i++) {
          if (vtx_mc_p == vtxMC_save[i]) {
            indexVtx = i;
            break;
  	  }
        }
        if (indexVtx >= 0) {
          printf("More than one MC-parent for MC-vertex %dn",indexVtx);
        }
        else {
          if (fVerbose > 1) printf("Create new T49VertexMCRootn");
          mcVertex = new T49VertexMCRoot();
          MCVertexList->Add(mcVertex);        
          mcVertex->SetX(vtx_mc_p->x);
          mcVertex->SetY(vtx_mc_p->y);
          mcVertex->SetZ(vtx_mc_p->z);
          if (fVerbose > 1) printf("Set parentn");
          mcVertex->SetParent(fNumMCTracks);
          mcTrack->SetStopVertex(fNumMCVertices);
          vtxMC_save[fNumMCVertices] = vtx_mc_p;
          fNumMCVertices++;
        }
      }
      else {
        if (fVerbose > 1) printf("No stop vertex foundn");
      }

      fNumMCTracks++;

    }

    if (fVerbose) {
      printf("Finished copying MC datan");
    }

  }

  if (!(fEventMC->GetNEvent() % 100) || fVerbose) {
    printf("Multiplicity: %d out of %d (%g)n"
          ,fNumTracks,TotalTracks,fEventMC->GetEveto());
    if (fFillSecondary) {
      printf("Secondary vertices (tracks): %d (%d)n"
            ,fNumSecVertices,fNumSecTracks);
    }
    if (fFillMCTracks) {
      printf("MC tracks: %d, MC vertices: %dn"
            ,fNumMCTracks,fNumMCVertices);
    }
  }

  delete [] rtrack_save;
  delete [] strack_save;
  delete [] vtx_save;
  delete [] vtxMC_save;

  return kTRUE;

}

//______________________________________________________________________________
 Bool_t TRootDSMC::FillMCEvent() 
{
  //
  //  Copies only the MC data for the current event.
  //

  Int_t          i;
  Int_t          num;
  Int_t          Dummy;   

  event_mc_t    *event_mc_p;
  track_mc_t    *track_mc_p;
  vertex_mc_t   *vtx_mc_p;

  T49ParticleMCRoot     *mcTrack;
  T49VertexMCRoot       *mcVertex;
  TClonesArray          *MCTrackArray;
  TObjArray             *MCVertexList;

  fNumMCTracks   = 0;
  fNumMCVertices = 0;

  if (fVerbose) {
    printf("Copy Monte Carlo datan");
  }

  // The TClonesArray for the MC tracks
  if (!(MCTrackArray  = fEventMC->GetMCParticles())) {
    printf("No Monte Carlo track array found in eventn");
  }
  // The TObjArray for the MC vertices
  if (!(MCVertexList  = fEventMC->GetMCVertices())) {
    printf("No Monte Carlo vertex list found in eventn");
  }

  event_mc_p = (event_mc_t *) get_objects("event_mc",&Dummy);
  if (event_mc_p == NULL) {
    cerr << "No event_mc_t in this event" << endl;
    return kFALSE;
  }

  // Used to create the links between the different lists
  Int_t         NVtxMCSave = event_mc_p->n_vertex + 1;
  vertex_mc_t **vtxMC_save = new vertex_mc_t*[NVtxMCSave];
  for (i = 0; i < NVtxMCSave;  i++) vtxMC_save[i] = NULL;

  if (fVerbose) {
    printf("Fill event_mc_t data (%d)n",fVerbose);
  }
  fEventMC->SetNRun(event_mc_p->n_run);
  fEventMC->SetNEvent(event_mc_p->n_event);

  // Get pointer to the track_mc_t structure 
  track_mc_p = (track_mc_t *) get_objects("track_mc",&num);

  // Loop over all MC tracks.
  for (Int_t inum = 0; inum < num; inum++) {

    // Create a new MC track and add it to the track array
    if (fVerbose > 1) printf("Create T49ParticleMCRootn");
    mcTrack = new (MCTrackArray->operator[](fNumMCTracks)) T49ParticleMCRoot();	

    // Fill the MC track
    FillMCTrack(mcTrack,track_mc_p);

    // Save the start-vertex
    if (fVerbose > 1) printf("Save the start vertexn");
    if ((vtx_mc_p = track_mc_p->start_vertex_p)) {
      Int_t indexVtx = -1;
      for (i = 0; i < fNumMCVertices; i++) {
        if (vtx_mc_p == vtxMC_save[i]) {
          indexVtx = i;
          break;
	}
      }
      if (indexVtx >= 0) {
        mcVertex = (T49VertexMCRoot *) MCVertexList->At(indexVtx);
        if (fVerbose > 1) printf("Vertex already copied %dn",indexVtx);
      }
      else {
       if (fNumMCVertices < NVtxMCSave) {   
         if (fVerbose > 1) printf("Create new T49VertexMCRootn");
          mcVertex = new T49VertexMCRoot();
          MCVertexList->Add(mcVertex);        
          mcVertex->SetX(vtx_mc_p->x);
          mcVertex->SetY(vtx_mc_p->y);
          mcVertex->SetZ(vtx_mc_p->z);
          vtxMC_save[fNumMCVertices] = vtx_mc_p;
          indexVtx = fNumMCVertices;
          fNumMCVertices++;
          if (fNumMCVertices == NVtxMCSave) {
            printf("More than %d MC vertices, skipping rest of verticesn",NVtxMCSave);
	  }
        }
        else {
          mcVertex = 0;
        }
      }
      if      (mcVertex->GetFirstDaughter()  < 0) {
        if (fVerbose > 1) printf("Set first daughtern");
        mcVertex->SetFirstDaughter(fNumMCTracks);
        mcTrack->SetStartVertex(indexVtx);
      }
      else if (mcVertex->GetSecondDaughter() < 0) {
        if (fVerbose > 1) printf("Set second daughtern");
        mcVertex->SetSecondDaughter(fNumMCTracks);
        mcTrack->SetStartVertex(indexVtx);
      }
      else if (mcVertex->GetThirdDaughter()  < 0) {
        if (fVerbose > 1) printf("Set third daughtern");
        mcVertex->SetThirdDaughter(fNumMCTracks);
        mcTrack->SetStartVertex(indexVtx);
      }
      else if (mcVertex->GetForthDaughter()  < 0) {
        if (fVerbose > 1) printf("Set forth daughtern");
        mcVertex->SetForthDaughter(fNumMCTracks);
        mcTrack->SetStartVertex(indexVtx);
      }
      else {
        if (fVerbose) printf("More than four MC-daughters for MC-vertex %dn"
                            ,indexVtx);
        mcTrack->SetStartVertex(indexVtx);
      }
    }
    else {
      if (fVerbose > 1) printf("No start vertex foundn");
    }

    // Save the stop vertex
    if (fVerbose > 1) printf("Save the stop vertexn");
    if ((vtx_mc_p = track_mc_p->stop_vertex_p)) {
      Int_t indexVtx = -1;
      for (i = 0; i < fNumMCVertices; i++) {
        if (vtx_mc_p == vtxMC_save[i]) {
          indexVtx = i;
          break;
        }
      }
      if (indexVtx >= 0) {
        printf("More than one MC-parent for MC-vertex %dn",indexVtx);
      }
      else {
        if (fNumMCVertices < NVtxMCSave) {   
          if (fVerbose > 1) printf("Create new T49VertexMCRootn");
          mcVertex = new T49VertexMCRoot();
          MCVertexList->Add(mcVertex);        
          mcVertex->SetX(vtx_mc_p->x);
          mcVertex->SetY(vtx_mc_p->y);
          mcVertex->SetZ(vtx_mc_p->z);
          if (fVerbose > 1) printf("Set parentn");
          mcVertex->SetParent(fNumMCTracks);
          mcTrack->SetStopVertex(fNumMCVertices);
          vtxMC_save[fNumMCVertices] = vtx_mc_p;
          fNumMCVertices++;
          if (fNumMCVertices == NVtxMCSave) {
            printf("More than %d MC vertices, skipping rest of verticesn",NVtxMCSave);
	  }
        }
      }
    }
    else {
      if (fVerbose > 1) printf("No stop vertex foundn");
    }

    track_mc_p++;
    fNumMCTracks++;

  }

  if (fVerbose) {
    printf("Finished copying MC datan");
    printf("MC tracks: %d, MC vertices: %dn"
          ,fNumMCTracks,fNumMCVertices);
  }

  delete [] vtxMC_save;

  return kTRUE;

}

//______________________________________________________________________________
 Bool_t TRootDSMC::FillMCTrack(T49ParticleMCRoot *mcTrack, track_mc_t *track_mc_p) 
{
  //
  // Copy data from a track_mc_t structure into a T49ParticleMCRoot object.
  //

  point_tof_mc_t *point_tof_mc_p;
  point_tpc_mc_t *point_mc;
  Float_t xfirst = 0;
  Float_t yfirst = 0;
  Float_t zfirst = 0;
  Float_t xlast  = 0;
  Float_t ylast  = 0;
  Float_t zlast  = 0;
 
  // The track_mc_t information
  mcTrack->SetPid(track_mc_p->pid);
  mcTrack->SetPx(track_mc_p->px);
  mcTrack->SetPy(track_mc_p->py);
  mcTrack->SetPz(track_mc_p->pz);
  mcTrack->SetNPoint(track_mc_p->n_vt1_hit,0);
  mcTrack->SetNPoint(track_mc_p->n_vt2_hit,1);
  mcTrack->SetNPoint(track_mc_p->n_mtl_hit
                   + track_mc_p->n_mtr_hit,2);
  mcTrack->SetNPoint(track_mc_p->n_vt1_hit
                   + track_mc_p->n_vt2_hit
                   + track_mc_p->n_mtl_hit
                   + track_mc_p->n_mtr_hit,3);

  // First and last MC points
  if (track_mc_p->pnt_vt1_p) {
    point_mc = track_mc_p->pnt_vt1_p;
    xlast  = point_mc->x;
    ylast  = point_mc->y;
    zlast  = point_mc->z;
    mcTrack->SetXLast(0,point_mc->x);
    mcTrack->SetYLast(0,point_mc->y);
    mcTrack->SetZLast(0,point_mc->z);
    while (point_mc->next_trk_pnt_p) {
      point_mc = point_mc->next_trk_pnt_p;
    }
    xfirst = point_mc->x;
    yfirst = point_mc->y;
    zfirst = point_mc->z;
    mcTrack->SetXFirst(0,point_mc->x);
    mcTrack->SetYFirst(0,point_mc->y);
    mcTrack->SetZFirst(0,point_mc->z);
  }
  else {
    mcTrack->SetXLast(0,0);
    mcTrack->SetYLast(0,0);
    mcTrack->SetZLast(0,0);
    mcTrack->SetXFirst(0,0);
    mcTrack->SetYFirst(0,0);
    mcTrack->SetZFirst(0,0);
  }
  if (track_mc_p->pnt_vt2_p) {
    point_mc = track_mc_p->pnt_vt2_p;
    xlast  = point_mc->x;
    ylast  = point_mc->y;
    zlast  = point_mc->z;
    mcTrack->SetXLast(1,point_mc->x);
    mcTrack->SetYLast(1,point_mc->y);
    mcTrack->SetZLast(1,point_mc->z);
    while (point_mc->next_trk_pnt_p) {
      point_mc = point_mc->next_trk_pnt_p;
    }
    if (xfirst == 0) {
      xfirst = point_mc->x;
      yfirst = point_mc->y;
      zfirst = point_mc->z;
    }
    mcTrack->SetXFirst(1,point_mc->x);
    mcTrack->SetYFirst(1,point_mc->y);
    mcTrack->SetZFirst(1,point_mc->z);
  }
  else {
    mcTrack->SetXLast(1,0);
    mcTrack->SetYLast(1,0);
    mcTrack->SetZLast(1,0);
    mcTrack->SetXFirst(1,0);
    mcTrack->SetYFirst(1,0);
    mcTrack->SetZFirst(1,0);
  }
  if ((track_mc_p->pnt_mtl_p) || 
      (track_mc_p->pnt_mtr_p)) {
    point_mc = track_mc_p->pnt_mtl_p 
             ? track_mc_p->pnt_mtl_p 
             : track_mc_p->pnt_mtr_p;
    xlast  = point_mc->x;
    ylast  = point_mc->y;
    zlast  = point_mc->z;
    mcTrack->SetXLast(2,point_mc->x);
    mcTrack->SetYLast(2,point_mc->y);
    mcTrack->SetZLast(2,point_mc->z);
    while (point_mc->next_trk_pnt_p) {
      point_mc=point_mc->next_trk_pnt_p;
    }
    if (xfirst == 0) {
      xfirst = point_mc->x;
      yfirst = point_mc->y;
      zfirst = point_mc->z;
    }
    mcTrack->SetXFirst(2,point_mc->x);
    mcTrack->SetYFirst(2,point_mc->y);
    mcTrack->SetZFirst(2,point_mc->z);
  }
  else {
    mcTrack->SetXLast(2,0);
    mcTrack->SetYLast(2,0);
    mcTrack->SetZLast(2,0);
    mcTrack->SetXFirst(2,0);
    mcTrack->SetYFirst(2,0);
    mcTrack->SetZFirst(2,0);
  }
  mcTrack->SetXFirst(xfirst);
  mcTrack->SetYFirst(yfirst);
  mcTrack->SetZFirst(zfirst);
  mcTrack->SetXLast(xlast);
  mcTrack->SetYLast(ylast);
  mcTrack->SetZLast(zlast);

  // The TOF information
  if ((point_tof_mc_p = track_mc_p->pnt_tof_p)) {
    mcTrack->SetTofId(point_tof_mc_p->id_det);
    mcTrack->SetTofX(point_tof_mc_p->x);
    mcTrack->SetTofY(point_tof_mc_p->y);
    mcTrack->SetTofPathl(point_tof_mc_p->path_length);
  }

  return kTRUE;

}





ROOT page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.