// $Header: /data/reve-cvs/reve/alicore/VSDCreator.cxx,v 1.1.1.1 2005/12/04 19:42:49 aljam Exp $

#include "VSDCreator.h"

#include <Reve/TTreeTools.h>

#include <AliITSLoader.h>
#include <AliTPCTrackHitsV2.h>
#include <AliPDG.h>
#include <AliHit.h>
#include <AliSimDigits.h>
#include <AliKalmanTrack.h>
#include <AliESD.h>
#include <AliESDV0MI.h>
#include <AliTPCclusterMI.h>
#include <AliTPCClustersRow.h>
#include <AliITS.h>
#include <AliITSclusterV2.h>
#include <AliTrackReference.h>
#include <AliESDkink.h>

#include <AliRun.h>
#include <AliTPCParam.h>

#include <TSystem.h>

using namespace Reve;
using namespace AliReve;

using namespace std;

//______________________________________________________________________
// VSDCreator
//

ClassImp(VSDCreator);

 void VSDCreator::init()
{
  mKineType = KT_Standard;
  mDataDir  = ".";
  
  mTPCHitRes = 2;
  mTRDHitRes = 2;

  pRunLoader = 0;

  // Particles not in ROOT's PDG database occuring in ALICE
  AliPDG::AddParticlesToPdgDataBase();
  {
    TDatabasePDG *pdgDB = TDatabasePDG::Instance();
    // const Int_t kspe=50000000;
    const Int_t kion=10000000;

    const Double_t kAu2Gev=0.9314943228;
    const Double_t khSlash = 1.0545726663e-27;
    const Double_t kErg2Gev = 1/1.6021773349e-3;
    const Double_t khShGev = khSlash*kErg2Gev;
    const Double_t kYear2Sec = 3600*24*365.25;

    pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
		       0,1,"Ion",kion+10020);
    pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
		       khShGev/(12.33*kYear2Sec),1,"Ion",kion+10030);
    pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
		       khShGev/(12.33*kYear2Sec),2,"Ion",kion+20040);
    pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
		       0,2,"Ion",kion+20030);
  }

  // AliKalmanTrack::SetConvConst(1); 
}

 VSDCreator::VSDCreator()
{
  init();
}

 VSDCreator::VSDCreator(const Text_t* name, const Text_t* title) :
  VSD(name, title)
{
  init();
}

/**************************************************************************/

 void VSDCreator::CreateVSD(const Text_t* data_dir, Int_t event,
			   const Text_t* vsd_file)
{
  static const Exc_t _eh("VSDCreator::CreateVSD ");

  mDataDir = data_dir;
  mEvent   = event;

  string galice_file (Form("%s/galice.root", mDataDir.Data()));
 
  // printf("Acces file to open runloader %s n", mDataDir.Data());

  if(gSystem->AccessPathName(galice_file.c_str(), kReadPermission)) {
    throw(_eh + "Can not read file '" + galice_file + "'.");
  }
  pRunLoader = AliRunLoader::Open(galice_file.c_str());
  if(pRunLoader == 0)
    throw(_eh + "AliRunLoader::Open failed.");

  pRunLoader->LoadgAlice();
  Int_t status = pRunLoader->GetEvent(mEvent);
  if(status)
    throw(_eh + Form("GetEvent(%d) failed, exit code %s.", mEvent, status));

  pRunLoader->LoadHeader();
  pRunLoader->LoadKinematics();
  pRunLoader->LoadTrackRefs();
  pRunLoader->LoadHits();

  // GledNS::PushFD();

  TFile* file = TFile::Open(vsd_file, "RECREATE", "ALICE VisualizationDataSummary");
  mDirectory = new TDirectory("Event0", "");

  CreateTrees();

  file->Write();  
  file->Close();
  delete file; 
  mDirectory =0;

  //GledNS::PopFD();

  // clean after the VSD data was sucessfuly 
  // written
  mTreeK      = 0;
  mTreeH      = 0;
  //mTreeTR     = 0;
  mTreeC      = 0;
  mTreeV0     = 0;
  mTreeKK     = 0;
  mTreeR      = 0;
  mTreeGI     = 0;

  pRunLoader->UnloadAll();
  delete pRunLoader;
  if(gAlice) {
    delete gAlice; gAlice = 0;
  }
  pRunLoader = 0;
}

 void VSDCreator::CreateTrees()
{
  static const Exc_t _eh("VSDCreator::CreateTrees ");

  if(mDirectory == 0)
    throw(_eh + "output directory not set.");

  try {
    ConvertKinematics();
  } catch(Exc_t& exc) { warn_caller(exc); }

  try {
    ConvertHits();
  } catch(Exc_t& exc) { warn_caller(exc); }

  try {
    ConvertClusters();
  } catch(Exc_t& exc) { warn_caller(exc); }

  try {
    ConvertRecTracks();
  } catch(Exc_t& exc) {
    warn_caller(exc + " Skipping V0 extraction.");
    goto end_esd_processing;
  }

  try {
    ConvertV0();
  } catch(Exc_t& exc) { warn_caller(exc); }

  try {
    ConvertKinks();
  } catch(Exc_t& exc) { warn_caller(exc); }

 end_esd_processing:

  try {
    ConvertGenInfo();
  } catch(Exc_t& exc) { warn_caller(exc); }
   
}

/**************************************************************************/
// Kinematics
/**************************************************************************/

 void VSDCreator::ConvertKinematics()
{
  static const Exc_t _eh("VSDCreator::ConvertKinematics ");

  if(mTreeK != 0) 
    throw (_eh + "kinematics already converted");

  TTree* treek = pRunLoader->TreeK();
  if(treek == 0) {
    warn_caller(_eh + "no kinematics.");
    return;
  }
  TParticle tp, *_tp = &tp;
  treek->SetBranchAddress("Particles", &_tp);

  mDirectory->cd();
  mTreeK = new TTree("Kinematics", "TParticles sorted by Label");
 
  Text_t* prim_selection = 0;
  switch(mKineType) {
  case KT_Standard:     prim_selection = "fMother[0] == -1"; break;
  case KT_ProtonProton: prim_selection = "fStatusCode > 0";  break;
  }
  TTreeQuery evl;
  Int_t nprimary = evl.Select(treek, prim_selection);

  Int_t nentries = (Int_t)treek->GetEntries();
  vector<MCTrack>  vmc(nentries);
  for (Int_t idx=0; idx<nentries; idx++) {
    // index to entry
    Int_t ent = (idx < nprimary) ?
      idx + nentries - nprimary :
      idx - nprimary;

    treek->GetEntry(ent);
    vmc[idx]       = tp;
    vmc[idx].label = idx;
  }

  // read track refrences 
  TTree* mTreeTR =  pRunLoader->TreeTR();

  if(mTreeTR == 0) {
    warn_caller(_eh + "no TrackRefs; some data will not be available.");
  } else {
    TClonesArray* RunArrayTR = 0;
    mTreeTR->SetBranchAddress("AliRun", &RunArrayTR);

    Int_t nPrimaries = (Int_t) mTreeTR->GetEntries();
    for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
      // printf("START mTreeTR->GetEntry(%d) n",iPrimPart);
      mTreeTR->GetEntry(iPrimPart);
      // printf("END mTreeTR->GetEntry(%d) n",iPrimPart);
    
      for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
	AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef); 
	Int_t track = trackRef->GetTrack();
	if(track < nentries && track > 0){ 
	  MCTrack& mct = vmc[track];	
	  if(trackRef->TestBit(kNotDeleted)) {
	    mct.decayed   = true;
	    mct.t_decay   = trackRef->GetTime();
	    mct.V_decay.x = trackRef->X();
	    mct.V_decay.y = trackRef->Y();
	    mct.V_decay.z = trackRef->Z();
	    mct.P_decay.x = trackRef->Px();
	    mct.P_decay.y = trackRef->Py();
	    mct.P_decay.z = trackRef->Pz();
	    if(TMath::Abs(mct.GetPdgCode()) == 11)
	      mct.decayed = false; // a bug in TreeTR
	  }
	}       
      }
    }
  }

  mTreeK->Branch("K", "Reve::MCTrack",  &mpK, fBuffSize);

  for(vector<MCTrack>::iterator k=vmc.begin(); k!=vmc.end(); ++k) {
    MCTrack& mct = *k;
    mK = mct;

    TParticle* m  = &mct;
    Int_t      mi = mct.label;
    while(m->GetMother(0) != -1) {
      mi = m->GetMother(0);
      m = &vmc[mi];
    }
    mK.eva_label = mi;

    mTreeK->Fill();
  }

  mTreeK->BuildIndex("label");
}

/**************************************************************************/
// Hits
/**************************************************************************/

namespace {

  struct Detector 
  {
    const char*   name;
    const char*   hitbranch;
    unsigned char detidx;
  };

  Detector detects[] = {
    { "ITS",  "AliITShit",         0 },
    { "TPC",  "AliTPCTrackHitsV2", 1 },
    { "TRD",  "AliTRDhit",         2 },
    { "TOF",  "AliTOFhit",         3 },
    // { "RICH", "AliRICHhit",        4 },
    { 0 }
  };

}

/**************************************************************************/

 void VSDCreator::ConvertHits()
{
  static const Exc_t _eh("VSDCreator::ConvertHits ");

  if(mTreeH != 0)
    throw(_eh + "hits already converted.");

  mDirectory->cd();
  mTreeH =  new TTree("Hits", "Combined detector hits.");
  mTreeH->Branch("H", "Reve::Hit", &mpH, fBuffSize);
 
  map<Int_t, Int_t> hmap;
  // parameters for ITS, TPC hits filtering
  Float_t x,y,z, x1,y1,z1;
  Float_t tpc_sqr_res = mTPCHitRes*mTPCHitRes;
  Float_t trd_sqr_res = mTRDHitRes*mTRDHitRes;

  int l=0;
  // load hits from the rest of detectors
  while(detects[l].name != 0) {
    Detector& det = detects[l++];

    switch(det.detidx) {
    case 1: { 
      Int_t count = 0;
      TTree* treeh = pRunLoader->GetTreeH(det.name, false);
      if(treeh == 0) {
	warn_caller(_eh + "no hits for "+ det.name +".");
	continue;
      }
      AliTPCTrackHitsV2 hv2, *_hv2=&hv2; 
      treeh->SetBranchAddress("TPC2", &_hv2);
      Int_t np = treeh->GetEntries();
      for(Int_t i=0; i<np; i++){
	treeh->GetEntry(i);
	Int_t eva_idx = np -i -1;
	if (hv2.First() == 0) continue;
        x = y = z = 0;
	do {
	  AliHit* ah = hv2.GetHit();
	  x1=ah->X();y1=ah->Y();z1=ah->Z();
	  if((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1) > tpc_sqr_res) {
	    mH.det_id    = det.detidx;
	    mH.subdet_id = 0;
	    mH.label     = ah->Track();
	    mH.eva_label = eva_idx;
	    mH.V.x = x1; mH.V.y = y1; mH.V.z = z1;
	    mTreeH->Fill();
	    hmap[mH.label]++;
	    x = x1; y = y1; z = z1;
	    count++;
	  }
	} while (hv2.Next());
      }
      // printf("%d entries in TPChits n",count);
      break;
    }
    default: {
      TTree* treeh = pRunLoader->GetTreeH(det.name, false);
      if(treeh == 0) {
	warn_caller(_eh + "no hits for "+ det.name +".");
	continue;
      }
      TClonesArray *arr = new TClonesArray(det.hitbranch);
      treeh->SetBranchAddress(det.name, &arr);
      Int_t np = treeh->GetEntries();
      // in TreeH files hits are grouped in clones arrays
      // each eva particle has its own clone array 
      for (Int_t i=0; i<np; i++) {
	treeh->GetEntry(i);
	Int_t eva_idx = np -i -1;
	Int_t nh=arr->GetEntriesFast();
	x = y = z = 0;
	// printf("%d entry %d hits for primary %d n", i, nh, eva_idx);
	for (Int_t j=0; j<nh; j++) {
	  AliHit* ali_hit = (AliHit*)arr->UncheckedAt(j);
	  mH.det_id    = det.detidx;
	  mH.subdet_id = 0;
	  mH.label     = ali_hit->GetTrack();
	  mH.eva_label = eva_idx;
	  mH.V.Set(ali_hit->X(), ali_hit->Y(), ali_hit->Z());
	  if(det.detidx == 2) {
	    x1=ali_hit->X();y1=ali_hit->Y();z1=ali_hit->Z();
	    if((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1) < trd_sqr_res) continue;
	    x=x1; y=y1; z=z1;
	  } 
	  hmap[mH.label]++;
	  mTreeH->Fill(); 
	}
      }
      delete arr;
      break;
    } // end default 
    } // end switch
  } // end while
  

  //set geninfo
  for(map<Int_t, Int_t>::iterator j=hmap.begin(); j!=hmap.end(); ++j) {
    get_geninfo(j->first)->n_hits += j->second;
  }
  
}

/**************************************************************************/
// Clusters
/**************************************************************************/

 void VSDCreator::ConvertClusters()
{
  static const Exc_t _eh("VSDCreator::ConvertClusters ");

  if(mTreeC != 0)
    throw(_eh + "clusters already converted.");

  mDirectory->cd();
  mTreeC =  new TTree("Clusters", "rec clusters");
  mTreeC->Branch("C", "Reve::Cluster", &mpC, fBuffSize);

  try {
    ConvertITSClusters();
  } catch(Exc_t& exc) { warn_caller(exc); }

  try {
    ConvertTPCClusters();
  } catch(Exc_t& exc) { warn_caller(exc); }
}

/**************************************************************************/

 void VSDCreator::ConvertTPCClusters()
{
  static const Exc_t _eh("VSDCreator::ConvertTPCClusters ");

  auto_ptr<TFile> f 
    ( TFile::Open(Form("%s/TPC.RecPoints.root", mDataDir.Data())) );
  if(!f.get())
    throw(_eh + "can not open 'TPC.RecPoints.root' file.");
    
  auto_ptr<TDirectory> d
    ( (TDirectory*) f->Get(Form("Event%d", mEvent)) );
  if(!d.get())
    throw(_eh + Form("event directory '%d' not found.", 0));

  auto_ptr<TTree> tree( (TTree*) d->Get("TreeR") );
  if(!tree.get())
    throw(_eh + "'TreeR' not found.");

  auto_ptr<AliTPCParam> par( get_tpc_param(_eh) );

  AliTPCClustersRow  clrow, *_clrow=&clrow;
  AliTPCclusterMI   *cl;
  _clrow->SetClass("AliTPCclusterMI");
  tree->SetBranchAddress("Segment", &_clrow);

  // count clusters
  Int_t nClusters = 0;
  Int_t n_ent = tree->GetEntries();
  for (Int_t n=0; n<n_ent; n++) {
    tree->GetEntry(n);
    nClusters += _clrow->GetArray()->GetEntriesFast();
  }

  // calculate xyz for a cluster and add it to container 
  Double_t x,y,z;
  Float_t cs, sn, tmp;
  map<Int_t, Int_t> cmap;

  for (Int_t n=0; n<tree->GetEntries(); n++) {
    tree->GetEntry(n);
    Int_t ncl = _clrow->GetArray()->GetEntriesFast();
    if(ncl > 0) {
      Int_t sec,row;
      par->AdjustSectorRow(_clrow->GetID(),sec,row);    
      while (ncl--) {
	if(_clrow->GetArray()) {
	  // cl = new AliTPCclusterMI(*(AliTPCclusterMI*)_clrow->GetArray()->UncheckedAt(ncl));
	  cl = (AliTPCclusterMI*)_clrow->GetArray()->UncheckedAt(ncl);
          if(cl->GetLabel(0) >= 0){
	    x = par->GetPadRowRadii(sec,row); y = cl->GetY(); z = cl->GetZ();
	    par->AdjustCosSin(sec,cs,sn);
	    tmp = x*cs-y*sn; y= x*sn+y*cs; x=tmp; 

	    mC.det_id    = 1;
	    mC.subdet_id = 0;
	    mC.label[0]  = cl->GetLabel(0);
	    mC.label[1]  = cl->GetLabel(1);
	    mC.label[2]  = cl->GetLabel(2);
	    mC.V.Set(x, y, z);

	    mTreeC->Fill();
	    { int i = 0;
	      while(i < 3 && mC.label[i])
		cmap[mC.label[i++]]++;
	    }
	  }
	}
      }
    }
  }
  //set geninfo
  for(map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j) {
    get_geninfo(j->first)->n_clus += j->second;
  }
}

/**************************************************************************/

 void VSDCreator::ConvertITSClusters()
{
  static const Exc_t _eh("VSDCreator::ConvertITSClusters ");

  auto_ptr<TFile> f 
    ( TFile::Open(Form("%s/ITS.RecPoints.root", mDataDir.Data())) );
  if(!f.get())
    throw(_eh + "can not open 'ITS.RecPoints.root' file.");
    
  auto_ptr<TDirectory> d
    ( (TDirectory*) f->Get(Form("Event%d", mEvent)) );
  if(!d.get())
    throw(_eh + Form("event directory '%d' not found.", 0));

  auto_ptr<TTree> tree( (TTree*) d->Get("TreeR") );
  if(!tree.get())
    throw(_eh + "'TreeR' not found.");

  AliITSLoader* ITSld =  (AliITSLoader*) pRunLoader->GetLoader("ITSLoader");
  //AliITS* pITS = ITSld->GetITS();
  AliITSgeom* geom = ITSld->GetITSgeom();
  //AliITSgeom* geom = new AliITSgeom();
  //geom->ReadNewFile("/home/aljam/ITSgeometry.det");

  //printf("alice ITS geom %p n",geom );

  if(!geom)
    throw(_eh + "can not find ITS geometry");

  TClonesArray *arr = new TClonesArray("AliITSclusterV2");
  tree->SetBranchAddress("Clusters", &arr);
  Int_t nmods = tree->GetEntries();
  Float_t gc[3];
  map<Int_t, Int_t> cmap;

  for (Int_t mod=0; mod<nmods; mod++) {
    tree->GetEntry(mod);
    Int_t nc=arr->GetEntriesFast();
    for (Int_t j=0; j<nc; j++) {
      AliITSclusterV2* recp = (AliITSclusterV2*)arr->UncheckedAt(j);

      Double_t rot[9];     
      geom->GetRotMatrix(mod,rot);
      Int_t lay,lad,det;   
      geom->GetModuleId(mod,lay,lad,det);
      Float_t tx,ty,tz;    
      geom->GetTrans(lay,lad,det,tx,ty,tz);     

      Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
      Double_t phi1=TMath::Pi()/2+alpha;
      if(lay==1) phi1+=TMath::Pi();

      Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
      Float_t  r=tx*cp+ty*sp;
      gc[0]= r*cp - recp->GetY()*sp;
      gc[1]= r*sp + recp->GetY()*cp;
      gc[2]= recp->GetZ();

      mC.det_id    = 0;
      mC.subdet_id = 0;
      mC.label[0]  = recp->GetLabel(0);
      mC.label[1]  = recp->GetLabel(1);
      mC.label[2]  = recp->GetLabel(2);
      mC.V.x       = r*cp - recp->GetY()*sp;
      mC.V.y       = r*sp + recp->GetY()*cp;
      mC.V.z       = recp->GetZ();
      mTreeC->Fill();
      { int i = 0;
	while(i < 3 && mC.label[i])
	  cmap[mC.label[i++]]++;
      }
    } 

    for(map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j) {
      get_geninfo(j->first)->n_clus += j->second;
    }
  }
  delete arr;
}

/**************************************************************************/
// ESD
/**************************************************************************/

 void VSDCreator::ConvertRecTracks()
{
  static const Exc_t _eh("VSDCreator::ConvertRecTracks ");

  if(mTreeR != 0)
    throw(_eh + "tracks already converted.");

  mDirectory->cd();
  mTreeR =  new TTree("RecTracks", "rec tracks");

  mTreeR->Branch("R", "Reve::RecTrack", &mpR, 512*1024,1);
 
  TFile f(Form("%s/AliESDs.root", mDataDir.Data()));
  if(!f.IsOpen())
    throw(_eh + "no AliESDs.root file.");

  TTree* tree = (TTree*) f.Get("esdTree");
  if (tree == 0) 
    throw(_eh + "no esdTree.");

 
  AliESD *fEvent=0;  
  tree->SetBranchAddress("ESD", &fEvent);
  tree->GetEntry(mEvent); 

 
  // reconstructed tracks
  AliESDtrack* esd_t;
  Double_t     dbuf[3];
  for (Int_t n=0; n<fEvent->GetNumberOfTracks(); n++) {
    esd_t = fEvent->GetTrack(n);

    mR.label  = esd_t->GetLabel();
    mR.status = (Int_t) esd_t->GetStatus();
    mR.sign   = esd_t->GetSign();
    esd_t->GetXYZ(dbuf);    mR.V.Set(dbuf);
    esd_t->GetPxPyPz(dbuf); mR.P.Set(dbuf);
    Double_t ep = esd_t->GetP();
    mR.beta = ep/TMath::Sqrt(ep*ep + TMath::C()*TMath::C()*esd_t->GetMass()*esd_t->GetMass());
    mTreeR->Fill();
  }
  mTreeR->BuildIndex("label");
}

/**************************************************************************/

 void VSDCreator::ConvertV0()
{
  static const Exc_t _eh("VSDCreator::ConvertV0 ");

  if(mTreeV0 != 0)
    throw(_eh + "V0 already converted.");

  mDirectory->cd();
  mTreeV0 =  new TTree("V0", "V0 points");

  mTreeV0->Branch("V0", "Reve::RecV0", &mpV0, 512*1024,1);

  TFile f(Form("%s/AliESDs.root", mDataDir.Data()));
  if(!f.IsOpen()){
    throw(_eh + "no AliESDs.root file.");
  }

  TTree* tree = (TTree*) f.Get("esdTree");
  if (tree == 0) 
    throw(_eh + "no esdTree.");

  AliESD *fEvent=0;  
  tree->SetBranchAddress("ESD", &fEvent);
  tree->GetEntry(mEvent); 

  for (Int_t n =0; n< fEvent->GetNumberOfV0MIs(); n++) {
    AliESDV0MI* av = fEvent->GetV0MI(n);
    Double_t x,y,cos,sin; 

    mV0.status = av->GetStatus();
    // Point of closest approach
    mV0.V_ca.x = av->GetXr(0); 
    mV0.V_ca.y = av->GetXr(1);
    mV0.V_ca.z = av->GetXr(2);
    // set birth vertex of neutral particle     
    { Double_t p[3];
      av->GetXYZ(p[0], p[1], p[2]);
      mV0.V0_birth.Set(p);
    }

    // momentum and position of negative particle
    mV0.P_neg.Set(av->GetPMp());
    // read AliExternalTrackParam
    x = av->GetParamM()->X(); 
    y = av->GetParamM()->Y(); 
    cos = TMath::Cos( av->GetParamM()->Alpha());
    sin = TMath::Sin( av->GetParamM()->Alpha());
    mV0.V_neg.x = x*cos - y*sin;
    mV0.V_neg.y = x*sin + y*cos;
    mV0.V_neg.z = av->GetParamM()->Z();

    // momentum and position of positive particle
    mV0.P_pos.Set(av->GetPPp());
    x = av->GetParamP()->X();
    y = av->GetParamP()->Y();
    cos = TMath::Cos(av->GetParamP()->Alpha());
    sin = TMath::Sin(av->GetParamP()->Alpha());
    mV0.V_pos.x = x*cos - y*sin;
    mV0.V_pos.y = x*sin + y*cos;
    mV0.V_pos.z = av->GetParamP()->Z();

    mV0.label = 0; // !!!! mother label unknown
    mV0.pdg   = av->GetPdgCode();
    // daughter indices
    mV0.d_label[0] = av->fLab[0];
    mV0.d_label[1] = av->fLab[1];

    // printf("V0 convert labels(%d,%d) index(%d,%d)n", 
    //	   av->GetLab(0), av->GetLab(1),
    //	   av->GetIndex(0), av->GetIndex(1));

    mTreeV0->Fill();
  }
  // if(fEvent->GetNumberOfV0MIs()) mTreeV0->BuildIndex("label");
}

/**************************************************************************/

 void VSDCreator::ConvertKinks()
{
  static const Exc_t _eh("VSDCreator::ConvertKinks ");

  if(mTreeKK != 0)
    throw(_eh + "Kinks already converted.");

  mDirectory->cd();
  mTreeKK =  new TTree("Kinks", "ESD Kinks");

  mTreeKK->Branch("KK", "Reve::RecKink", &mpKK, fBuffSize);

  TFile f(Form("%s/AliESDs.root", mDataDir.Data()));
  if(!f.IsOpen()){
    throw(_eh + "no AliESDs.root file.");
  }

  TTree* tree = (TTree*) f.Get("esdTree");
  if (tree == 0) 
    throw(_eh + "no esdTree.");

  AliESD *fEvent=0;  
  tree->SetBranchAddress("ESD", &fEvent);
  tree->GetEntry(mEvent); 

  //  printf("CONVERT KINK Read %d entries in tree kinks n",  fEvent->GetNumberOfKinks());
  for (Int_t n =0; n< fEvent->GetNumberOfKinks(); n++) {
    AliESDkink* kk = fEvent->GetKink(n);

    Double_t x,y,cos,sin; 

    mKK.label  = kk->fLab[0];
    mKK.status = Int_t(kk->fStatus);

    // reconstructed kink position
    mKK.label_sec = kk->fLab[1];
    mKK.V_kink.Set(&kk->fXr[0]);

    // momentum and position of mother 
    mKK.P.x = kk->fParamMother.Px();
    mKK.P.y = kk->fParamMother.Py();
    mKK.P.z = kk->fParamMother.Pz();
    const Double_t* par =  kk->fParamMother.GetParameter();
    // printf("KINK Pt %f, %f n",1/kk->fParamMother.Pt(),par[4] );
    mKK.sign = (par[4] < 0) ? -1 : 1;

    x = kk->fParamMother.X(); 
    y = kk->fParamMother.Y(); 
    cos = TMath::Cos( kk->fParamMother.Alpha());
    sin = TMath::Sin( kk->fParamMother.Alpha());
    mKK.V.z = x*cos - y*sin;
    mKK.V.y = x*sin + y*cos;
    mKK.V.z = kk->fParamMother.Z();
   
    // momentum and position of daughter 
    mKK.P_sec.x = kk->fParamDaughter.Px();
    mKK.P_sec.y = kk->fParamDaughter.Py();
    mKK.P_sec.z = kk->fParamDaughter.Pz();

    x = kk->fParamDaughter.X(); 
    y = kk->fParamDaughter.Y(); 
    cos = TMath::Cos(kk->fParamDaughter.Alpha());
    sin = TMath::Sin( kk->fParamDaughter.Alpha());
    mKK.V_end.x = x*cos - y*sin;
    mKK.V_end.y = x*sin + y*cos;
    mKK.V_end.z = kk->fParamDaughter.Z();

    mTreeKK->Fill();
  }
  if(fEvent->GetNumberOfKinks()) mTreeKK->BuildIndex("label");
}
/**************************************************************************/
// GenInfo
/**************************************************************************/

 void VSDCreator::ConvertGenInfo()
{
  static const Exc_t _eh("VSDCreator::ConvertGenInfo ");

  if(mTreeGI != 0)
    throw(_eh + "GI already converted.");

  mDirectory->cd();
  mTreeGI = new TTree("GenInfo", "Objects prepared for cross querry");

  GenInfo::Class()->IgnoreTObjectStreamer(true);
  mTreeGI->Branch("GI", "Reve::GenInfo",  &mpGI, fBuffSize);
  mTreeGI->Branch("K.", "Reve::MCTrack",  &mpK);
  mTreeGI->Branch("R.", "Reve::RecTrack", &mpR);

  for(map<Int_t, GenInfo*>::iterator j=mGenInfoMap.begin(); j!=mGenInfoMap.end(); ++j) {
    mGI = *(j->second);
    mGI.label = j->first;
    mTreeK->GetEntry(j->first);

    if(mTreeR) {
      Int_t re = mTreeR->GetEntryNumberWithIndex(j->first);
      if(re != -1) 
	mGI.is_rec = true;
    }
    //    Int_t has_v0 =  mTreeV0->GetEntryNumberWithIndex(j->first);
    //if (has_v0 != -1)
    //  mGI.has_V0 = true;
    if (mTreeKK) {
      Int_t has_kk =  mTreeKK->GetEntryNumberWithIndex(j->first);
      if (has_kk != -1)
	mGI.has_kink = true;
    }
    mTreeGI->Fill();
  }
  mGenInfoMap.clear();
}

/**************************************************************************/
/**************************************************************************/
// Protected methods
/**************************************************************************/
/**************************************************************************/

 AliTPCParam* VSDCreator::get_tpc_param(const Exc_t& eh)
{
  auto_ptr<TFile> fp( TFile::Open(Form("%s/galice.root", mDataDir.Data())) );
  if(!fp.get())
    throw(eh + "can not open 'galice.root' file.");
  AliTPCParam* par = (AliTPCParam *) fp->Get("75x40_100x60_150x60");
  if(!par)
    throw(eh + "TPC data not found.");
  return par;
}



 GenInfo* VSDCreator::get_geninfo(Int_t label)
{
  // printf("get_geninfo %dn", label);
  GenInfo* gi;
  map<Int_t, GenInfo*>::iterator i = mGenInfoMap.find(label);
  if(i == mGenInfoMap.end()) {
    gi =  new GenInfo();
    mGenInfoMap[label] = gi;
  } else {
    gi = i->second;
  }
  return gi;
}