[u/mrichter/AliRoot.git] / RAW / AliHoughFilter.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// high level filter algorithm for TPC using a hough transformation          //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TStopwatch.h>

#include "AliHLTStandardIncludes.h"
#include "AliHLTLogging.h"
#include "AliHLTTransform.h"
#include "AliHLTHough.h"
#include "AliLog.h"
#include <AliHLTITSclusterer.h>
#include <AliHLTITSVertexerZ.h>
#include <AliHLTITStracker.h>

#include "AliHoughFilter.h"

#include "AliRawReaderRoot.h"
#include <AliMagF.h>
#include <AliKalmanTrack.h>
#include <AliITSgeom.h>
#include <AliESDVertex.h>

ClassImp(AliHoughFilter)

//_____________________________________________________________________________
AliHoughFilter::AliHoughFilter():
fPtmin(0),
fITSgeom(NULL)
{
// default constructor

  // Init debug level
  AliHLTLog::fgLevel = AliHLTLog::kError;
  if (AliDebugLevel() > 0) AliHLTLog::fgLevel = AliHLTLog::kWarning;
  if (AliDebugLevel() > 1) AliHLTLog::fgLevel = AliHLTLog::kInformational;
  if (AliDebugLevel() > 2) AliHLTLog::fgLevel = AliHLTLog::kDebug;

  // Init TPC HLT geometry
  const char *path = gSystem->Getenv("ALICE_ROOT");
  Char_t pathname[1024];
  strcpy(pathname,path);
  strcat(pathname,"/HLT/src");
  if (!AliHLTTransform::Init(pathname, kFALSE))
    AliError("HLT initialization failed!");

  // Init magnetic field
  AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance();
  AliTracker::SetFieldMap(field,kTRUE);
  fPtmin = 0.1*AliHLTTransform::GetSolenoidField();

  // Init ITS geometry
  fITSgeom = new AliITSgeom();
  fITSgeom->ReadNewFile("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymmFMD.det");
  if (!fITSgeom)
    AliError("ITS geometry not found!");

  // Init PID
  AliPID pid;
}

//_____________________________________________________________________________
Bool_t AliHoughFilter::Filter(AliRawVEvent* event, AliESDEvent* esd)
{
  // Run fast online reconstruction
  // based on the HLT tracking algorithms

  TStopwatch globaltimer;
  globaltimer.Start();

  TStopwatch timer;
  timer.Start();

  TTree *treeITSclusters = new TTree("TreeL3ITSclusters"," "); //make a tree
  treeITSclusters->SetDirectory(0);

  RunITSclusterer(event,treeITSclusters);
  RunITSvertexer(esd,treeITSclusters);
  RunTPCtracking(event,esd);
  RunITStracking(esd,treeITSclusters);
  delete treeITSclusters;

  AliInfo(Form("Event filter has finished in %f seconds\n\n\n\n\n\n",globaltimer.RealTime()));

  return kTRUE;
}

//_____________________________________________________________________________
void AliHoughFilter::RunITSclusterer(AliRawVEvent* event, TTree *treeClusters)
{
  // Run ITS Clusterer
  // The clusters are stored in a tree

  TStopwatch timer;
  timer.Start();

  AliHLTITSclusterer clusterer(0);
  AliRawReader *itsrawreader=new AliRawReaderRoot(event);
  clusterer.Digits2Clusters(itsrawreader,treeClusters);
  delete itsrawreader;
  AliInfo(Form("ITS clusterer has finished in %f seconds\n",timer.RealTime()));

}


//_____________________________________________________________________________
void AliHoughFilter::RunITSvertexer(AliESDEvent* esd, TTree *treeClusters)
{
  // Run SPD vertexerZ
  // Store the result in the ESD

  TStopwatch timer;
  timer.Start();

  AliHLTITSVertexerZ vertexer;
  AliESDVertex *vertex = vertexer.FindVertexForCurrentEvent(fITSgeom,treeClusters);
  esd->SetVertex(vertex);
  AliInfo(Form("ITS vertexer has finished in %f seconds\n",timer.RealTime()));

}

//_____________________________________________________________________________
void AliHoughFilter::RunTPCtracking(AliRawVEvent* event, AliESDEvent* esd)
{
  // Run hough transform tracking in TPC
  // The z of the vertex is taken from the ESD
  // The result of the tracking is stored in the ESD
  TStopwatch timer;
  timer.Start();

  const AliESDVertex *vertex = esd->GetVertex();
  Float_t zvertex = vertex->GetZv();

  AliHLTHough *hough1 = new AliHLTHough();
    
  hough1->SetThreshold(4);
  hough1->CalcTransformerParams(fPtmin);
  hough1->SetPeakThreshold(70,-1);
  hough1->Init(100,4,event,zvertex);
  hough1->SetAddHistograms();

  AliHLTHough *hough2 = new AliHLTHough();
  
  hough2->SetThreshold(4);
  hough2->CalcTransformerParams(fPtmin);
  hough2->SetPeakThreshold(70,-1);
  hough2->Init(100,4,event,zvertex);
  hough2->SetAddHistograms();

  Int_t nglobaltracks = 0;
  /* In case we run HLT code in 2 threads */
  hough1->StartProcessInThread(0,17);
  hough2->StartProcessInThread(18,35);

  if(hough1->WaitForThreadFinish())
    ::Fatal("AliHLTHough::WaitForThreadFinish"," Can not join the required thread! ");
  if(hough2->WaitForThreadFinish())
    ::Fatal("AliHLTHough::WaitForThreadFinish"," Can not join the required thread! ");

    /* In case we run HLT code in the main thread
    for(Int_t slice=0; slice<=17; slice++)
      {
	hough1->ReadData(slice,0);
	hough1->Transform();
	hough1->AddAllHistogramsRows();
	hough1->FindTrackCandidatesRow();
	hough1->AddTracks();
      }
    for(Int_t slice=18; slice<=35; slice++)
      {
	hough2->ReadData(slice,0);
	hough2->Transform();
	hough2->AddAllHistogramsRows();
	hough2->FindTrackCandidatesRow();
	hough2->AddTracks();
      }
    */

  nglobaltracks += hough1->FillESD(esd);
  nglobaltracks += hough2->FillESD(esd);

    /* In case we want to debug the ESD
    gSystem->MakeDirectory("hough1");
    hough1->WriteTracks("./hough1");
    gSystem->MakeDirectory("hough2");
    hough2->WriteTracks("./hough2");
    */

  delete hough1;
  delete hough2;

  AliInfo(Form("\nHough Transformer has found %d TPC tracks in %f seconds\n\n", nglobaltracks,timer.RealTime()));

}

//_____________________________________________________________________________
void AliHoughFilter::RunITStracking(AliESDEvent* esd, TTree *treeClusters)
{
  // Run the ITS tracker
  // The tracks from the HT TPC tracking are used as seeds
  // The prologated tracks are updated in the ESD
  TStopwatch timer;
  timer.Start();

  Double_t vtxPos[3];
  Double_t vtxErr[3]={0.005,0.005,0.010};
  const AliESDVertex *vertex = esd->GetVertex();
  vertex->GetXYZ(vtxPos);

  AliHLTITStracker itsTracker(0);
  itsTracker.SetVertex(vtxPos,vtxErr);

  itsTracker.LoadClusters(treeClusters);
  itsTracker.Clusters2Tracks(esd);
  itsTracker.UnloadClusters();
  AliInfo(Form("ITS tracker has finished in %f seconds\n",timer.RealTime()));

}
Commit	Line	Data
	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// high level filter algorithm for TPC using a hough transformation //
	21	// //
	22	///////////////////////////////////////////////////////////////////////////////
	23
	24	#include <TStopwatch.h>
	25
	26	#include "AliHLTStandardIncludes.h"
	27	#include "AliHLTLogging.h"
	28	#include "AliHLTTransform.h"
	29	#include "AliHLTHough.h"
	30	#include "AliLog.h"
	31	#include <AliHLTITSclusterer.h>
	32	#include <AliHLTITSVertexerZ.h>
	33	#include <AliHLTITStracker.h>
	34
	35	#include "AliHoughFilter.h"
	36
	37	#include "AliRawReaderRoot.h"
	38	#include <AliMagF.h>
	39	#include <AliKalmanTrack.h>
	40	#include <AliITSgeom.h>
	41	#include <AliESDVertex.h>
	42
	43	ClassImp(AliHoughFilter)
	44
	45	//_____________________________________________________________________________
	46	AliHoughFilter::AliHoughFilter():
	47	fPtmin(0),
	48	fITSgeom(NULL)
	49	{
	50	// default constructor
	51
	52	// Init debug level
	53	AliHLTLog::fgLevel = AliHLTLog::kError;
	54	if (AliDebugLevel() > 0) AliHLTLog::fgLevel = AliHLTLog::kWarning;
	55	if (AliDebugLevel() > 1) AliHLTLog::fgLevel = AliHLTLog::kInformational;
	56	if (AliDebugLevel() > 2) AliHLTLog::fgLevel = AliHLTLog::kDebug;
	57
	58	// Init TPC HLT geometry
	59	const char *path = gSystem->Getenv("ALICE_ROOT");
	60	Char_t pathname[1024];
	61	strcpy(pathname,path);
	62	strcat(pathname,"/HLT/src");
	63	if (!AliHLTTransform::Init(pathname, kFALSE))
	64	AliError("HLT initialization failed!");
	65
	66	// Init magnetic field
	67	AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance();
	68	AliTracker::SetFieldMap(field,kTRUE);
	69	fPtmin = 0.1*AliHLTTransform::GetSolenoidField();
	70
	71	// Init ITS geometry
	72	fITSgeom = new AliITSgeom();
	73	fITSgeom->ReadNewFile("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymmFMD.det");
	74	if (!fITSgeom)
	75	AliError("ITS geometry not found!");
	76
	77	// Init PID
	78	AliPID pid;
	79	}
	80
	81	//_____________________________________________________________________________
	82	Bool_t AliHoughFilter::Filter(AliRawVEvent* event, AliESDEvent* esd)
	83	{
	84	// Run fast online reconstruction
	85	// based on the HLT tracking algorithms
	86
	87	TStopwatch globaltimer;
	88	globaltimer.Start();
	89
	90	TStopwatch timer;
	91	timer.Start();
	92
	93	TTree *treeITSclusters = new TTree("TreeL3ITSclusters"," "); //make a tree
	94	treeITSclusters->SetDirectory(0);
	95
	96	RunITSclusterer(event,treeITSclusters);
	97	RunITSvertexer(esd,treeITSclusters);
	98	RunTPCtracking(event,esd);
	99	RunITStracking(esd,treeITSclusters);
	100	delete treeITSclusters;
	101
	102	AliInfo(Form("Event filter has finished in %f seconds\n\n\n\n\n\n",globaltimer.RealTime()));
	103
	104	return kTRUE;
	105	}
	106
	107	//_____________________________________________________________________________
	108	void AliHoughFilter::RunITSclusterer(AliRawVEvent* event, TTree *treeClusters)
	109	{
	110	// Run ITS Clusterer
	111	// The clusters are stored in a tree
	112
	113	TStopwatch timer;
	114	timer.Start();
	115
	116	AliHLTITSclusterer clusterer(0);
	117	AliRawReader *itsrawreader=new AliRawReaderRoot(event);
	118	clusterer.Digits2Clusters(itsrawreader,treeClusters);
	119	delete itsrawreader;
	120	AliInfo(Form("ITS clusterer has finished in %f seconds\n",timer.RealTime()));
	121
	122	}
	123
	124
	125	//_____________________________________________________________________________
	126	void AliHoughFilter::RunITSvertexer(AliESDEvent* esd, TTree *treeClusters)
	127	{
	128	// Run SPD vertexerZ
	129	// Store the result in the ESD
	130
	131	TStopwatch timer;
	132	timer.Start();
	133
	134	AliHLTITSVertexerZ vertexer;
	135	AliESDVertex *vertex = vertexer.FindVertexForCurrentEvent(fITSgeom,treeClusters);
	136	esd->SetVertex(vertex);
	137	AliInfo(Form("ITS vertexer has finished in %f seconds\n",timer.RealTime()));
	138
	139	}
	140
	141	//_____________________________________________________________________________
	142	void AliHoughFilter::RunTPCtracking(AliRawVEvent* event, AliESDEvent* esd)
	143	{
	144	// Run hough transform tracking in TPC
	145	// The z of the vertex is taken from the ESD
	146	// The result of the tracking is stored in the ESD
	147	TStopwatch timer;
	148	timer.Start();
	149
	150	const AliESDVertex *vertex = esd->GetVertex();
	151	Float_t zvertex = vertex->GetZv();
	152
	153	AliHLTHough *hough1 = new AliHLTHough();
	154
	155	hough1->SetThreshold(4);
	156	hough1->CalcTransformerParams(fPtmin);
	157	hough1->SetPeakThreshold(70,-1);
	158	hough1->Init(100,4,event,zvertex);
	159	hough1->SetAddHistograms();
	160
	161	AliHLTHough *hough2 = new AliHLTHough();
	162
	163	hough2->SetThreshold(4);
	164	hough2->CalcTransformerParams(fPtmin);
	165	hough2->SetPeakThreshold(70,-1);
	166	hough2->Init(100,4,event,zvertex);
	167	hough2->SetAddHistograms();
	168
	169	Int_t nglobaltracks = 0;
	170	/* In case we run HLT code in 2 threads */
	171	hough1->StartProcessInThread(0,17);
	172	hough2->StartProcessInThread(18,35);
	173
	174	if(hough1->WaitForThreadFinish())
	175	::Fatal("AliHLTHough::WaitForThreadFinish"," Can not join the required thread! ");
	176	if(hough2->WaitForThreadFinish())
	177	::Fatal("AliHLTHough::WaitForThreadFinish"," Can not join the required thread! ");
	178
	179	/* In case we run HLT code in the main thread
	180	for(Int_t slice=0; slice<=17; slice++)
	181	{
	182	hough1->ReadData(slice,0);
	183	hough1->Transform();
	184	hough1->AddAllHistogramsRows();
	185	hough1->FindTrackCandidatesRow();
	186	hough1->AddTracks();
	187	}
	188	for(Int_t slice=18; slice<=35; slice++)
	189	{
	190	hough2->ReadData(slice,0);
	191	hough2->Transform();
	192	hough2->AddAllHistogramsRows();
	193	hough2->FindTrackCandidatesRow();
	194	hough2->AddTracks();
	195	}
	196	*/
	197
	198	nglobaltracks += hough1->FillESD(esd);
	199	nglobaltracks += hough2->FillESD(esd);
	200
	201	/* In case we want to debug the ESD
	202	gSystem->MakeDirectory("hough1");
	203	hough1->WriteTracks("./hough1");
	204	gSystem->MakeDirectory("hough2");
	205	hough2->WriteTracks("./hough2");
	206	*/
	207
	208	delete hough1;
	209	delete hough2;
	210
	211	AliInfo(Form("\nHough Transformer has found %d TPC tracks in %f seconds\n\n", nglobaltracks,timer.RealTime()));
	212
	213	}
	214
	215	//_____________________________________________________________________________
	216	void AliHoughFilter::RunITStracking(AliESDEvent* esd, TTree *treeClusters)
	217	{
	218	// Run the ITS tracker
	219	// The tracks from the HT TPC tracking are used as seeds
	220	// The prologated tracks are updated in the ESD
	221	TStopwatch timer;
	222	timer.Start();
	223
	224	Double_t vtxPos[3];
	225	Double_t vtxErr[3]={0.005,0.005,0.010};
	226	const AliESDVertex *vertex = esd->GetVertex();
	227	vertex->GetXYZ(vtxPos);
	228
	229	AliHLTITStracker itsTracker(0);
	230	itsTracker.SetVertex(vtxPos,vtxErr);
	231
	232	itsTracker.LoadClusters(treeClusters);
	233	itsTracker.Clusters2Tracks(esd);
	234	itsTracker.UnloadClusters();
	235	AliInfo(Form("ITS tracker has finished in %f seconds\n",timer.RealTime()));
	236
	237	}