[u/mrichter/AliRoot.git] / EVE / alice-macros / esd_tracks.C

// $Id$
// Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007

/**************************************************************************
 * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
 * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for          *
 * full copyright notice.                                                 *
 **************************************************************************/

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <TString.h>
#include <TSystem.h>
#include <TROOT.h>
#include <TMath.h>
#include <TGListTree.h>

#include <TEveVSDStructs.h>
#include <TEveManager.h>
#include <TEveTrackPropagator.h>

#include <EveBase/AliEveTrack.h>
#include <EveBase/AliEveTrackCounter.h>
#include <EveBase/AliEveMagField.h>
#include <EveBase/AliEveEventManager.h>

#include <AliESDEvent.h>
#include <AliESDtrackCuts.h>
#include <AliESDtrack.h>
#include <AliExternalTrackParam.h>

#include <AliPWG0Helper.h>

#endif


// Use inner-tpc track params when its refit failed.
Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;

// Use magnetic-field as retrieved from GRP.
Bool_t g_esd_tracks_true_field = kTRUE;

// Use Runge-Kutta track stepper.
Bool_t g_esd_tracks_rk_stepper = kFALSE;


//==============================================================================

void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
				Float_t magF, Float_t maxR)
{
  if (g_esd_tracks_true_field)
  {
    trkProp->SetMagFieldObj(new AliEveMagField);
  }
  else
  {
    trkProp->SetMagField(magF);
  }
  if (g_esd_tracks_rk_stepper)
  {
    trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
  }
  trkProp->SetMaxR(maxR);
}

//==============================================================================

TString esd_track_title(AliESDtrack* t)
{
  TString s;

  Int_t label = t->GetLabel(), index = t->GetID();
  TString idx(index == kMinInt ? "<undef>" : Form("%d", index));
  TString lbl(label == kMinInt ? "<undef>" : Form("%d", label));

  Double_t p[3], v[3];
  t->GetXYZ(v);
  t->GetPxPyPz(p);
  Double_t pt    = t->Pt();
  Double_t ptsig = TMath::Sqrt(t->GetSigma1Pt2());
  Double_t ptsq  = pt*pt;
  Double_t ptm   = pt / (1.0 + pt*ptsig);
  Double_t ptM   = pt / (1.0 - pt*ptsig);

  s = Form("Index=%s, Label=%s\nChg=%d, Pdg=%d\n"
	   "pT = %.3f + %.3f - %.3f [%.3f]\n"
           "P  = (%.3f, %.3f, %.3f)\n"
           "V  = (%.3f, %.3f, %.3f)\n",
	   idx.Data(), lbl.Data(), t->Charge(), 0,
	   pt, ptM - pt, pt - ptm, ptsig*ptsq,
           p[0], p[1], p[2],
           v[0], v[1], v[2]);

  Int_t   o;
  s += "Det (in,out,refit,pid):\n";
  o  = AliESDtrack::kITSin;
  s += Form("ITS (%d,%d,%d,%d)  ",  t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  o  = AliESDtrack::kTPCin;
  s += Form("TPC(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  o  = AliESDtrack::kTRDin;
  s += Form("TRD(%d,%d,%d,%d) ",    t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  o  = AliESDtrack::kTOFin;
  s += Form("TOF(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  o  = AliESDtrack::kHMPIDout;
  s += Form("HMPID(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
  s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));

  if (t->IsOn(AliESDtrack::kESDpid))
  {
    Double_t pid[5];
    t->GetESDpid(pid);
    s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
  }

  return s;
}

//==============================================================================

void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
{
  // Add additional track parameters as a path-mark to track.

  if (tp == 0)
    return;

  Double_t pbuf[3], vbuf[3];
  tp->GetXYZ(vbuf);
  tp->GetPxPyPz(pbuf);

  TEvePathMark pm(TEvePathMark::kReference);
  pm.fV.Set(vbuf);
  pm.fP.Set(pbuf);
  track->AddPathMark(pm);
}

//==============================================================================

AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
{
  // Make a standard track representation and put it into given container.

  // Choose which parameters to use a track's starting point.
  // If gkFixFailedITSExtr is TRUE (FALSE by default) and
  // if ITS refit failed, take track parameters at inner TPC radius.

  const AliExternalTrackParam* tp = at;

  Bool_t innerTaken = kFALSE;
  if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
  {
    tp = at->GetInnerParam();
    innerTaken = kTRUE;
  }

  AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
  track->SetAttLineAttMarker(cont);
  track->SetName(Form("AliEveTrack %d", at->GetID()));
  track->SetElementTitle(esd_track_title(at));
  track->SetSourceObject(at);

  // Add inner/outer track parameters as path-marks.
  if (at->IsOn(AliESDtrack::kTPCrefit))
  {
    if ( ! innerTaken)
    {
      esd_track_add_param(track, at->GetInnerParam());
    }
    esd_track_add_param(track, at->GetOuterParam());
  }

  return track;
}


//==============================================================================
// esd_tracks()
//==============================================================================

TEveTrackList* esd_tracks()
{
  AliESDEvent* esd = AliEveEventManager::AssertESD();

  TEveTrackList* cont = new TEveTrackList("ESD Tracks");
  cont->SetMainColor(6);

  esd_track_propagator_setup(cont->GetPropagator(),
			     0.1*esd->GetMagneticField(), 520);

  gEve->AddElement(cont);

  Int_t count = 0;
  for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
  {
    ++count;
    AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);

    cont->AddElement(track);
  }
  cont->SetTitle(Form("N=%d", count));
  cont->MakeTracks();

  gEve->Redraw3D();

  return cont;
}

TEveTrackList* esd_tracks_MI()
{
  AliESDEvent* esd = AliEveEventManager::AssertESD();

  TEveTrackList* cont = new TEveTrackList("ESD Tracks MI");
  cont->SetLineColor(5);
  gEve->AddElement(cont);

  Int_t count = 0;
  for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
  {
    ++count;
    AliESDtrack* at = esd->GetTrack(n);
    AliEveTrack* l = new AliEveTrack(at, cont->GetPropagator());
    l->SetName(Form("ESDTrackMI %d", at->GetID()));
    l->SetElementTitle(esd_track_title(at));
    l->SetAttLineAttMarker(cont);
    l->SetSourceObject(at);

    at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);

    l->SetLockPoints(kTRUE);
    cont->AddElement(l);
  }
  cont->SetTitle(Form("N=%d", count));

  gEve->Redraw3D();

  return cont;
}


//==============================================================================
// esd_tracks_from_array()
//==============================================================================

TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
{
  // Retrieves AliESDTrack's from collection.
  // See example usage with AliAnalysisTrackCuts in the next function.

  if (esd == 0) esd = AliEveEventManager::AssertESD();

  TEveTrackList* cont = new TEveTrackList("ESD Tracks");
  cont->SetMainColor(6);

  esd_track_propagator_setup(cont->GetPropagator(),
			     0.1*esd->GetMagneticField(), 520);

  gEve->AddElement(cont);

  Int_t    count = 0;
  TIter    next(col);
  TObject *obj;
  while ((obj = next()) != 0)
  {
    if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
    {
      Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
	      obj->GetName(), obj->GetTitle());
      continue;
    }

    ++count;
    AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);

    AliEveTrack* track = esd_make_track(at, cont);
    cont->AddElement(track);
  }
  cont->SetTitle(Form("N=%d", count));
  cont->MakeTracks();

  gEve->Redraw3D();

  return cont;
}

void esd_tracks_alianalcuts_demo()
{
  AliESDEvent* esd = AliEveEventManager::AssertESD();

  AliESDtrackCuts atc;
  atc.SetPtRange(0.1, 5);
  atc.SetRapRange(-1, 1);

  esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
}


//==============================================================================
// esd_tracks_by_category
//==============================================================================

Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
{
  // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
  // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)

  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];
  esdTrack->GetImpactParameters(b,bCov);
  if (bCov[0] <= 0 || bCov[2] <= 0)
  {
    printf("Estimated b resolution lower or equal zero!\n");
    bCov[0] = bCov[2] = 0;
  }
  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);

  // -----------------------------------
  // How to get to a n-sigma cut?
  //
  // The accumulated statistics from 0 to d is
  //
  // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
  // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
  //
  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
  // Can this be expressed in a different way?

  if (bRes[0] == 0 || bRes[1] == 0)
    return -1;

  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));

  // stupid rounding problem screws up everything:
  // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
  if (TMath::Exp(-d * d / 2) < 1e-10)
    return 1000;

  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  return d;
}

TEveElementList* esd_tracks_by_category()
{
   // Import ESD tracks, separate them into several containers
   // according to primary-vertex cut and ITS&TPC refit status.

   AliESDEvent* esd = AliEveEventManager::AssertESD();

   TEveElementList* cont = new TEveElementList("ESD Tracks by category");
   gEve->AddElement(cont);

   const Int_t   nCont = 9;
   const Float_t maxR  = 520;
   const Float_t magF  = 0.1*esd->GetMagneticField();

   TEveTrackList *tl[nCont];
   Int_t          tc[nCont];
   Int_t          count = 0;

   tl[0] = new TEveTrackList("Sigma < 3");
   tc[0] = 0;
   esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
   tl[0]->SetMainColor(3);
   cont->AddElement(tl[0]);

   tl[1] = new TEveTrackList("3 < Sigma < 5");
   tc[1] = 0;
   esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
   tl[1]->SetMainColor(7);
   cont->AddElement(tl[1]);

   tl[2] = new TEveTrackList("5 < Sigma");
   tc[2] = 0;
   esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
   tl[2]->SetMainColor(46);
   cont->AddElement(tl[2]);

   tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
   tc[3] = 0;
   esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
   tl[3]->SetMainColor(41);
   cont->AddElement(tl[3]);

   tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
   tc[4] = 0;
   esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
   tl[4]->SetMainColor(48);
   cont->AddElement(tl[4]);

   tl[5] = new TEveTrackList("no TPC refit");
   tc[5] = 0;
   esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
   tl[5]->SetMainColor(kRed);
   cont->AddElement(tl[5]);

   tl[6] = new TEveTrackList("ITS ncl>=3 & SPD Inner");
   tc[6] = 0;
   esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
   tl[6]->SetMainColor(kGreen);
   cont->AddElement(tl[6]);

   tl[7] = new TEveTrackList("ITS ncl>=3 & b<3 cm");
   tc[7] = 0;
   esd_track_propagator_setup(tl[7]->GetPropagator(), magF, maxR);
   tl[7]->SetMainColor(kMagenta - 9);
   cont->AddElement(tl[7]);

   tl[8] = new TEveTrackList("ITS others");
   tc[8] = 0;
   esd_track_propagator_setup(tl[8]->GetPropagator(), magF, maxR);
   tl[8]->SetMainColor(kRed + 2);
   cont->AddElement(tl[8]);

   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
   {
      AliESDtrack* at = esd->GetTrack(n);

      Float_t s  = get_sigma_to_vertex(at);
      Int_t   ti;
      if      (s <  3) ti = 0;
      else if (s <= 5) ti = 1;
      else             ti = 2;

      Int_t    nclits;
      Double_t dtobeam;

      if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
      {
         UChar_t itsclmap = at->GetITSClusterMap();
         Bool_t  spdinner = (itsclmap & 3) != 0;

         nclits = 0;
         for (Int_t iter = 0; iter < 6; ++iter)
            if (itsclmap & (1 << iter)) nclits ++;

         Double_t xyz[3];
         at->GetXYZ(xyz);
         dtobeam = TMath::Hypot(xyz[0], xyz[1]);

         if ((nclits >= 3) && (spdinner))
            ti = 6;
         else if ((nclits >= 3) && (dtobeam < 3.0))
            ti = 7;
         else
            ti = 8;
      }
      else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
      {
         ti = 5;
      }
      else if ( ! at->IsOn(AliESDtrack::kITSrefit))
      {
         ti = (ti == 2) ? 4 : 3;
      }

      TEveTrackList* tlist = tl[ti];
      ++tc[ti];
      ++count;

      AliEveTrack* track = esd_make_track(at, tlist);
      if (ti<6)
         track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
      else
         track->SetName(Form("ESD Track idx=%d, dxy=%5.3f cl=%i", at->GetID(), dtobeam, nclits));
      tlist->AddElement(track);
   }

   for (Int_t ti = 0; ti < nCont; ++ti)
   {
      TEveTrackList* tlist = tl[ti];
      tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
      tlist->SetTitle(Form("N tracks=%d", tc[ti]));

      tlist->MakeTracks();

      //    Bool_t good_cont = ti <= 1;
      Bool_t good_cont = ((ti == 6) || (ti == 7));
      if (AliEveTrackCounter::IsActive())
      {
         AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
      }
      else
      {
         if ( ! good_cont)
            tlist->SetLineStyle(6);
      }
   }
   cont->SetTitle(Form("N all tracks = %d", count));
   // ??? The following does not always work:
   cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);

   gEve->Redraw3D();

   return cont;
}


//==============================================================================
// esd_tracks_by_anal_cuts
//==============================================================================

AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;

TEveElementList* esd_tracks_by_anal_cuts()
{
  AliESDEvent* esd = AliEveEventManager::AssertESD();

  if (g_esd_tracks_anal_cuts == 0)
  {
    gSystem->Load("libPWG0base");
    gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/CreateStandardCuts.C");
    Int_t mode = AliPWG0Helper::kTPC;
    if (TMath::Abs(esd->GetMagneticField()) > 0.01)
      mode |= AliPWG0Helper::kFieldOn;
    gROOT->ProcessLine(Form("g_esd_tracks_anal_cuts = CreateTrackCuts(%d, kFALSE)", mode));
  }

  TEveElementList* cont = new TEveElementList("ESD Tracks by Analysis Cuts");
  gEve->AddElement(cont);

  const Int_t   nCont = 2;
  const Float_t maxR  = 520;
  const Float_t magF  = 0.1*esd->GetMagneticField();

  TEveTrackList *tl[nCont];
  Int_t          tc[nCont];
  Int_t          count = 0;

  tl[0] = new TEveTrackList("Passed");
  tc[0] = 0;
  esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
  tl[0]->SetMainColor(3);
  cont->AddElement(tl[0]);

  tl[1] = new TEveTrackList("Rejected");
  tc[1] = 0;
  esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
  tl[1]->SetMainColor(kRed);
  cont->AddElement(tl[1]);

  for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
  {
    AliESDtrack* at = esd->GetTrack(n);

    Float_t s  = get_sigma_to_vertex(at);
    Int_t   ti = (g_esd_tracks_anal_cuts->AcceptTrack(at)) ? 0 : 1;

    TEveTrackList* tlist = tl[ti];
    ++tc[ti];
    ++count;

    AliEveTrack* track = esd_make_track(at, tlist);
    track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
    tlist->AddElement(track);
  }

  for (Int_t ti = 0; ti < nCont; ++ti)
  {
    TEveTrackList* tlist = tl[ti];
    tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
    tlist->SetTitle(Form("N tracks=%d", tc[ti]));

    tlist->MakeTracks();

    Bool_t good_cont = ti < 1;
    if (AliEveTrackCounter::IsActive())
    {
      AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
    }
    else
    {
      if ( ! good_cont)
        tlist->SetLineStyle(6);
    }
  }
  cont->SetTitle(Form("N all tracks = %d", count));
  // ??? The following does not always work:
  cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);

  gEve->Redraw3D();

  return cont;
}
Commit	Line	Data
5a5a1232	1	// $Id$
d810d0de	2	// Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
	3
	4	/**************************************************************************
	5	* Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
	6	* See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for *
51346b82	7	* full copyright notice. *
d810d0de	8	**************************************************************************/
5a5a1232	9
e9ece5b4	10	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
e9ece5b4	11
15c12442	12	#include <TString.h>
	13	#include <TSystem.h>
	14	#include <TROOT.h>
	15	#include <TMath.h>
	16	#include <TGListTree.h>
	17
	18	#include <TEveVSDStructs.h>
	19	#include <TEveManager.h>
	20	#include <TEveTrackPropagator.h>
	21
	22	#include <EveBase/AliEveTrack.h>
	23	#include <EveBase/AliEveTrackCounter.h>
	24	#include <EveBase/AliEveMagField.h>
	25	#include <EveBase/AliEveEventManager.h>
	26
	27	#include <AliESDEvent.h>
	28	#include <AliESDtrackCuts.h>
	29	#include <AliESDtrack.h>
	30	#include <AliExternalTrackParam.h>
	31
b264afb8	32	#include <AliPWG0Helper.h>
e9ece5b4	33
	34	#endif
	35
	36
29207369	37	// Use inner-tpc track params when its refit failed.
	38	Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
	39
6336e5d7	40	// Use magnetic-field as retrieved from GRP.
	41	Bool_t g_esd_tracks_true_field = kTRUE;
	42
	43	// Use Runge-Kutta track stepper.
	44	Bool_t g_esd_tracks_rk_stepper = kFALSE;
	45
	46
	47	//==============================================================================
	48
	49	void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
	50	Float_t magF, Float_t maxR)
	51	{
	52	if (g_esd_tracks_true_field)
	53	{
	54	trkProp->SetMagFieldObj(new AliEveMagField);
	55	}
	56	else
	57	{
	58	trkProp->SetMagField(magF);
	59	}
	60	if (g_esd_tracks_rk_stepper)
	61	{
	62	trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
	63	}
	64	trkProp->SetMaxR(maxR);
	65	}
	66
	67	//==============================================================================
	68
29207369	69	TString esd_track_title(AliESDtrack* t)
ad5abc55	70	{
ad5abc55	71	TString s;
29207369	72
	73	Int_t label = t->GetLabel(), index = t->GetID();
	74	TString idx(index == kMinInt ? "<undef>" : Form("%d", index));
	75	TString lbl(label == kMinInt ? "<undef>" : Form("%d", label));
	76
	77	Double_t p[3], v[3];
	78	t->GetXYZ(v);
	79	t->GetPxPyPz(p);
d6962cc4	80	Double_t pt = t->Pt();
	81	Double_t ptsig = TMath::Sqrt(t->GetSigma1Pt2());
	82	Double_t ptsq = pt*pt;
	83	Double_t ptm = pt / (1.0 + pt*ptsig);
	84	Double_t ptM = pt / (1.0 - pt*ptsig);
29207369	85
29207369	86	s = Form("Index=%s, Label=%s\nChg=%d, Pdg=%d\n"
d6962cc4	87	"pT = %.3f + %.3f - %.3f [%.3f]\n"
	88	"P = (%.3f, %.3f, %.3f)\n"
	89	"V = (%.3f, %.3f, %.3f)\n",
29207369	90	idx.Data(), lbl.Data(), t->Charge(), 0,
d6962cc4	91	pt, ptM - pt, pt - ptm, ptsig*ptsq,
	92	p[0], p[1], p[2],
	93	v[0], v[1], v[2]);
29207369	94
ad5abc55	95	Int_t o;
d6962cc4	96	s += "Det (in,out,refit,pid):\n";
ad5abc55	97	o = AliESDtrack::kITSin;
6b7fd3a4	98	s += Form("ITS (%d,%d,%d,%d) ", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55	99	o = AliESDtrack::kTPCin;
6b7fd3a4	100	s += Form("TPC(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55	101	o = AliESDtrack::kTRDin;
6b7fd3a4	102	s += Form("TRD(%d,%d,%d,%d) ", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55	103	o = AliESDtrack::kTOFin;
6b7fd3a4	104	s += Form("TOF(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55	105	o = AliESDtrack::kHMPIDout;
6b7fd3a4	106	s += Form("HMPID(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
ad5abc55	107	s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));
86f12f78	108
29207369	109	if (t->IsOn(AliESDtrack::kESDpid))
	110	{
	111	Double_t pid[5];
	112	t->GetESDpid(pid);
	113	s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
	114	}
51346b82	115
29207369	116	return s;
86f12f78	117	}
86f12f78	118
6336e5d7	119	//==============================================================================
6336e5d7	120
e9ece5b4	121	void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
c222eb1d	122	{
29207369	123	// Add additional track parameters as a path-mark to track.
29207369	124
c222eb1d	125	if (tp == 0)
	126	return;
	127
29207369	128	Double_t pbuf[3], vbuf[3];
c222eb1d	129	tp->GetXYZ(vbuf);
	130	tp->GetPxPyPz(pbuf);
	131
a15e6d7d	132	TEvePathMark pm(TEvePathMark::kReference);
	133	pm.fV.Set(vbuf);
	134	pm.fP.Set(pbuf);
c222eb1d	135	track->AddPathMark(pm);
	136	}
	137
6336e5d7	138	//==============================================================================
6336e5d7	139
e9ece5b4	140	AliEveTrack* esd_make_track(AliESDtrack at, TEveTrackList cont)
29207369	141	{
	142	// Make a standard track representation and put it into given container.
	143
	144	// Choose which parameters to use a track's starting point.
	145	// If gkFixFailedITSExtr is TRUE (FALSE by default) and
	146	// if ITS refit failed, take track parameters at inner TPC radius.
e9ece5b4	147
e9ece5b4	148	const AliExternalTrackParam* tp = at;
29207369	149
	150	Bool_t innerTaken = kFALSE;
	151	if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
	152	{
	153	tp = at->GetInnerParam();
	154	innerTaken = kTRUE;
	155	}
	156
50ac85c6	157	AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
29207369	158	track->SetAttLineAttMarker(cont);
0e33c639	159	track->SetName(Form("AliEveTrack %d", at->GetID()));
29207369	160	track->SetElementTitle(esd_track_title(at));
9b6216c0	161	track->SetSourceObject(at);
29207369	162
	163	// Add inner/outer track parameters as path-marks.
	164	if (at->IsOn(AliESDtrack::kTPCrefit))
	165	{
	166	if ( ! innerTaken)
	167	{
	168	esd_track_add_param(track, at->GetInnerParam());
	169	}
	170	esd_track_add_param(track, at->GetOuterParam());
	171	}
a15e6d7d	172
29207369	173	return track;
29207369	174	}
86f12f78	175
6336e5d7	176
	177	//==============================================================================
	178	// esd_tracks()
	179	//==============================================================================
	180
29207369	181	TEveTrackList* esd_tracks()
5a5a1232	182	{
d810d0de	183	AliESDEvent* esd = AliEveEventManager::AssertESD();
5a5a1232	184
51346b82	185	TEveTrackList* cont = new TEveTrackList("ESD Tracks");
fbc350a3	186	cont->SetMainColor(6);
6336e5d7	187
	188	esd_track_propagator_setup(cont->GetPropagator(),
	189	0.1*esd->GetMagneticField(), 520);
5a5a1232	190
84aff7a4	191	gEve->AddElement(cont);
5a5a1232	192
29207369	193	Int_t count = 0;
29207369	194	for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
d31ac39d	195	{
29207369	196	++count;
0e33c639	197	AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
5a5a1232	198
29207369	199	cont->AddElement(track);
	200	}
	201	cont->SetTitle(Form("N=%d", count));
	202	cont->MakeTracks();
5a5a1232	203
29207369	204	gEve->Redraw3D();
5a5a1232	205
29207369	206	return cont;
29207369	207	}
5a5a1232	208
06219f85	209	TEveTrackList* esd_tracks_MI()
29207369	210	{
29207369	211	AliESDEvent* esd = AliEveEventManager::AssertESD();
c222eb1d	212
06219f85	213	TEveTrackList* cont = new TEveTrackList("ESD Tracks MI");
06219f85	214	cont->SetLineColor(5);
29207369	215	gEve->AddElement(cont);
c222eb1d	216
29207369	217	Int_t count = 0;
	218	for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
	219	{
	220	++count;
	221	AliESDtrack* at = esd->GetTrack(n);
06219f85	222	AliEveTrack* l = new AliEveTrack(at, cont->GetPropagator());
29207369	223	l->SetName(Form("ESDTrackMI %d", at->GetID()));
29207369	224	l->SetElementTitle(esd_track_title(at));
06219f85	225	l->SetAttLineAttMarker(cont);
9b6216c0	226	l->SetSourceObject(at);
c222eb1d	227
06219f85	228	at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
	229
	230	l->SetLockPoints(kTRUE);
29207369	231	cont->AddElement(l);
5a5a1232	232	}
29207369	233	cont->SetTitle(Form("N=%d", count));
32e219c2	234
84aff7a4	235	gEve->Redraw3D();
5a5a1232	236
	237	return cont;
	238	}
	239
6336e5d7	240
6336e5d7	241	//==============================================================================
a8600b56	242	// esd_tracks_from_array()
6336e5d7	243	//==============================================================================
a8600b56	244
84aff7a4	245	TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
5a5a1232	246	{
	247	// Retrieves AliESDTrack's from collection.
	248	// See example usage with AliAnalysisTrackCuts in the next function.
	249
d810d0de	250	if (esd == 0) esd = AliEveEventManager::AssertESD();
86f12f78	251
51346b82	252	TEveTrackList* cont = new TEveTrackList("ESD Tracks");
fbc350a3	253	cont->SetMainColor(6);
6336e5d7	254
	255	esd_track_propagator_setup(cont->GetPropagator(),
	256	0.1*esd->GetMagneticField(), 520);
5a5a1232	257
84aff7a4	258	gEve->AddElement(cont);
5a5a1232	259
5a5a1232	260	Int_t count = 0;
	261	TIter next(col);
	262	TObject *obj;
84aff7a4	263	while ((obj = next()) != 0)
d31ac39d	264	{
29207369	265	if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
	266	{
	267	Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
5a5a1232	268	obj->GetName(), obj->GetTitle());
	269	continue;
	270	}
	271
	272	++count;
e9ece5b4	273	AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
5a5a1232	274
0e33c639	275	AliEveTrack* track = esd_make_track(at, cont);
29207369	276	cont->AddElement(track);
5a5a1232	277	}
6b7fd3a4	278	cont->SetTitle(Form("N=%d", count));
5a5a1232	279	cont->MakeTracks();
32e219c2	280
84aff7a4	281	gEve->Redraw3D();
5a5a1232	282
	283	return cont;
	284	}
	285
	286	void esd_tracks_alianalcuts_demo()
	287	{
d810d0de	288	AliESDEvent* esd = AliEveEventManager::AssertESD();
5a5a1232	289
fcfc8b01	290	AliESDtrackCuts atc;
5a5a1232	291	atc.SetPtRange(0.1, 5);
	292	atc.SetRapRange(-1, 1);
	293
fcfc8b01	294	esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
5a5a1232	295	}
a8600b56	296
6336e5d7	297
6336e5d7	298	//==============================================================================
29207369	299	// esd_tracks_by_category
6336e5d7	300	//==============================================================================
a8600b56	301
	302	Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
	303	{
	304	// Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
	305	// Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
	306
	307	Float_t b[2];
	308	Float_t bRes[2];
	309	Float_t bCov[3];
	310	esdTrack->GetImpactParameters(b,bCov);
29207369	311	if (bCov[0] <= 0 \|\| bCov[2] <= 0)
29207369	312	{
d31ac39d	313	printf("Estimated b resolution lower or equal zero!\n");
29207369	314	bCov[0] = bCov[2] = 0;
a8600b56	315	}
	316	bRes[0] = TMath::Sqrt(bCov[0]);
	317	bRes[1] = TMath::Sqrt(bCov[2]);
	318
	319	// -----------------------------------
	320	// How to get to a n-sigma cut?
	321	//
	322	// The accumulated statistics from 0 to d is
	323	//
	324	// -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
	325	// -> 1 - Exp(-d*2) for a 2-dim gauss (dd = dxdx + dydy != n_sigma)
	326	//
	327	// It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)ErfInv(1 - Exp((-x*2)/2)
	328	// Can this be expressed in a different way?
	329
29207369	330	if (bRes[0] == 0 \|\| bRes[1] == 0)
a8600b56	331	return -1;
	332
	333	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	334
	335	// stupid rounding problem screws up everything:
	336	// if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
	337	if (TMath::Exp(-d * d / 2) < 1e-10)
	338	return 1000;
	339
	340	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	341	return d;
	342	}
	343
29207369	344	TEveElementList* esd_tracks_by_category()
a8600b56	345	{
440a8bc6	346	// Import ESD tracks, separate them into several containers
	347	// according to primary-vertex cut and ITS&TPC refit status.
	348
	349	AliESDEvent* esd = AliEveEventManager::AssertESD();
	350
	351	TEveElementList* cont = new TEveElementList("ESD Tracks by category");
	352	gEve->AddElement(cont);
	353
	354	const Int_t nCont = 9;
	355	const Float_t maxR = 520;
	356	const Float_t magF = 0.1*esd->GetMagneticField();
	357
	358	TEveTrackList *tl[nCont];
	359	Int_t tc[nCont];
	360	Int_t count = 0;
	361
	362	tl[0] = new TEveTrackList("Sigma < 3");
	363	tc[0] = 0;
	364	esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
	365	tl[0]->SetMainColor(3);
	366	cont->AddElement(tl[0]);
	367
	368	tl[1] = new TEveTrackList("3 < Sigma < 5");
	369	tc[1] = 0;
	370	esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
	371	tl[1]->SetMainColor(7);
	372	cont->AddElement(tl[1]);
	373
	374	tl[2] = new TEveTrackList("5 < Sigma");
	375	tc[2] = 0;
	376	esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
	377	tl[2]->SetMainColor(46);
	378	cont->AddElement(tl[2]);
	379
	380	tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
	381	tc[3] = 0;
	382	esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
	383	tl[3]->SetMainColor(41);
	384	cont->AddElement(tl[3]);
	385
	386	tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
	387	tc[4] = 0;
	388	esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
	389	tl[4]->SetMainColor(48);
	390	cont->AddElement(tl[4]);
	391
	392	tl[5] = new TEveTrackList("no TPC refit");
	393	tc[5] = 0;
	394	esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
	395	tl[5]->SetMainColor(kRed);
	396	cont->AddElement(tl[5]);
	397
	398	tl[6] = new TEveTrackList("ITS ncl>=3 & SPD Inner");
	399	tc[6] = 0;
	400	esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
	401	tl[6]->SetMainColor(kGreen);
	402	cont->AddElement(tl[6]);
	403
	404	tl[7] = new TEveTrackList("ITS ncl>=3 & b<3 cm");
	405	tc[7] = 0;
	406	esd_track_propagator_setup(tl[7]->GetPropagator(), magF, maxR);
	407	tl[7]->SetMainColor(kMagenta - 9);
	408	cont->AddElement(tl[7]);
	409
410	tl[8] = new TEveTrackList("ITS others");
411	tc[8] = 0;
412	esd_track_propagator_setup(tl[8]->GetPropagator(), magF, maxR);
413	tl[8]->SetMainColor(kRed + 2);
414	cont->AddElement(tl[8]);
415
416	for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
417	{
418	AliESDtrack* at = esd->GetTrack(n);
419
420	Float_t s = get_sigma_to_vertex(at);
421	Int_t ti;
422	if (s < 3) ti = 0;
423	else if (s <= 5) ti = 1;
424	else ti = 2;
425
426	Int_t nclits;
427	Double_t dtobeam;
428
429	if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
430	{
431	UChar_t itsclmap = at->GetITSClusterMap();
432	Bool_t spdinner = (itsclmap & 3) != 0;
433
434	nclits = 0;
435	for (Int_t iter = 0; iter < 6; ++iter)
436	if (itsclmap & (1 << iter)) nclits ++;
437
438	Double_t xyz[3];
439	at->GetXYZ(xyz);
440	dtobeam = TMath::Hypot(xyz[0], xyz[1]);
441
442	if ((nclits >= 3) && (spdinner))
443	ti = 6;
444	else if ((nclits >= 3) && (dtobeam < 3.0))
445	ti = 7;
446	else
447	ti = 8;
448	}
449	else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
450	{
451	ti = 5;
452	}
453	else if ( ! at->IsOn(AliESDtrack::kITSrefit))
454	{
455	ti = (ti == 2) ? 4 : 3;
456	}
457
458	TEveTrackList* tlist = tl[ti];
459	++tc[ti];
460	++count;
461
462	AliEveTrack* track = esd_make_track(at, tlist);
463	if (ti<6)
464	track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
465	else
466	track->SetName(Form("ESD Track idx=%d, dxy=%5.3f cl=%i", at->GetID(), dtobeam, nclits));
467	tlist->AddElement(track);
468	}
469
470	for (Int_t ti = 0; ti < nCont; ++ti)
471	{
472	TEveTrackList* tlist = tl[ti];
473	tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
474	tlist->SetTitle(Form("N tracks=%d", tc[ti]));
475
476	tlist->MakeTracks();
477
478	// Bool_t good_cont = ti <= 1;
479	Bool_t good_cont = ((ti == 6) \|\| (ti == 7));
480	if (AliEveTrackCounter::IsActive())
481	{
482	AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
483	}
484	else
485	{
486	if ( ! good_cont)
487	tlist->SetLineStyle(6);
488	}
489	}
490	cont->SetTitle(Form("N all tracks = %d", count));
491	// ??? The following does not always work:
492	cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
493
494	gEve->Redraw3D();
495
496	return cont;
a8600b56	497	}
5320a363	498
	499
	500	//==============================================================================
	501	// esd_tracks_by_anal_cuts
	502	//==============================================================================
	503
	504	AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;
	505
	506	TEveElementList* esd_tracks_by_anal_cuts()
	507	{
	508	AliESDEvent* esd = AliEveEventManager::AssertESD();
	509
	510	if (g_esd_tracks_anal_cuts == 0)
	511	{
	512	gSystem->Load("libPWG0base");
	513	gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/CreateStandardCuts.C");
5065a65b	514	Int_t mode = AliPWG0Helper::kTPC;
5320a363	515	if (TMath::Abs(esd->GetMagneticField()) > 0.01)
5320a363	516	mode \|= AliPWG0Helper::kFieldOn;
15c12442	517	gROOT->ProcessLine(Form("g_esd_tracks_anal_cuts = CreateTrackCuts(%d, kFALSE)", mode));
5320a363	518	}
	519
	520	TEveElementList* cont = new TEveElementList("ESD Tracks by Analysis Cuts");
	521	gEve->AddElement(cont);
	522
	523	const Int_t nCont = 2;
	524	const Float_t maxR = 520;
	525	const Float_t magF = 0.1*esd->GetMagneticField();
	526
	527	TEveTrackList *tl[nCont];
	528	Int_t tc[nCont];
	529	Int_t count = 0;
	530
	531	tl[0] = new TEveTrackList("Passed");
	532	tc[0] = 0;
	533	esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
	534	tl[0]->SetMainColor(3);
	535	cont->AddElement(tl[0]);
	536
	537	tl[1] = new TEveTrackList("Rejected");
	538	tc[1] = 0;
	539	esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
	540	tl[1]->SetMainColor(kRed);
	541	cont->AddElement(tl[1]);
	542
	543	for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
	544	{
	545	AliESDtrack* at = esd->GetTrack(n);
	546
	547	Float_t s = get_sigma_to_vertex(at);
	548	Int_t ti = (g_esd_tracks_anal_cuts->AcceptTrack(at)) ? 0 : 1;
	549
	550	TEveTrackList* tlist = tl[ti];
	551	++tc[ti];
	552	++count;
	553
	554	AliEveTrack* track = esd_make_track(at, tlist);
	555	track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
	556	tlist->AddElement(track);
	557	}
	558
	559	for (Int_t ti = 0; ti < nCont; ++ti)
	560	{
	561	TEveTrackList* tlist = tl[ti];
	562	tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
	563	tlist->SetTitle(Form("N tracks=%d", tc[ti]));
	564
	565	tlist->MakeTracks();
	566
	567	Bool_t good_cont = ti < 1;
	568	if (AliEveTrackCounter::IsActive())
	569	{
	570	AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
	571	}
	572	else
	573	{
	574	if ( ! good_cont)
	575	tlist->SetLineStyle(6);
	576	}
	577	}
	578	cont->SetTitle(Form("N all tracks = %d", count));
	579	// ??? The following does not always work:
	580	cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
	581
582	gEve->Redraw3D();
583
584	return cont;
585	}