[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 "AliL3StandardIncludes.h"
#include "AliL3Logging.h"
#include "AliL3Transform.h"
#include "AliL3Hough.h"
#include "AliLog.h"
#include <AliL3ITSclusterer.h>
#include <AliL3ITSVertexerZ.h>
#include <AliL3ITStracker.h>

#include "AliHoughFilter.h"

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

ClassImp(AliHoughFilter)

//_____________________________________________________________________________
AliHoughFilter::AliHoughFilter()
{
// default constructor

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

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

  // Init magnetic field
  AliMagF* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
  AliTracker::SetFieldMap(field,kTRUE);
  fPtmin = 0.1*AliL3Transform::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(AliRawEvent* event, AliESD* 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(AliRawEvent* event, TTree *treeClusters)
{
  // Run ITS Clusterer
  // The clusters are stored in a tree

  TStopwatch timer;
  timer.Start();

  if(!fITSgeom)
    AliError("ITS geometry not created!");
  AliL3ITSclusterer clusterer(fITSgeom);
  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(AliESD* esd, TTree *treeClusters)
{
  // Run SPD vertexerZ
  // Store the result in the ESD

  TStopwatch timer;
  timer.Start();

  AliL3ITSVertexerZ 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(AliRawEvent* event, AliESD* 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();

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

  AliL3Hough *hough2 = new AliL3Hough();
  
  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("AliL3Hough::WaitForThreadFinish"," Can not join the required thread! ");
  if(hough2->WaitForThreadFinish())
    ::Fatal("AliL3Hough::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(AliESD* 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);

  AliL3ITStracker itsTracker(fITSgeom);
  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
e10815f1	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
e10815f1	24	#include <TStopwatch.h>
	25
	26	#include "AliL3StandardIncludes.h"
	27	#include "AliL3Logging.h"
	28	#include "AliL3Transform.h"
	29	#include "AliL3Hough.h"
	30	#include "AliLog.h"
d6186083	31	#include <AliL3ITSclusterer.h>
	32	#include <AliL3ITSVertexerZ.h>
	33	#include <AliL3ITStracker.h>
e10815f1	34
	35	#include "AliHoughFilter.h"
	36
d6186083	37	#include "AliRawReaderRoot.h"
	38	#include <AliMagF.h>
	39	#include <AliMagFMaps.h>
	40	#include <AliKalmanTrack.h>
	41	#include <AliITSgeom.h>
e10815f1	42
	43	ClassImp(AliHoughFilter)
	44
	45	//_____________________________________________________________________________
	46	AliHoughFilter::AliHoughFilter()
	47	{
	48	// default constructor
	49
d6186083	50	// Init debug level
e10815f1	51	AliL3Log::fgLevel = AliL3Log::kError;
	52	if (AliDebugLevel() > 0) AliL3Log::fgLevel = AliL3Log::kWarning;
	53	if (AliDebugLevel() > 1) AliL3Log::fgLevel = AliL3Log::kInformational;
	54	if (AliDebugLevel() > 2) AliL3Log::fgLevel = AliL3Log::kDebug;
d6186083	55
	56	// Init TPC HLT geometry
	57	const char *path = gSystem->Getenv("ALICE_ROOT");
	58	Char_t pathname[1024];
	59	strcpy(pathname,path);
	60	strcat(pathname,"/HLT/src");
	61	if (!AliL3Transform::Init(pathname, kFALSE))
e10815f1	62	AliError("HLT initialization failed!");
d6186083	63
d6186083	64	// Init magnetic field
d6186083	65	AliMagF* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
c84a5e9e	66	AliTracker::SetFieldMap(field,kTRUE);
d6186083	67	fPtmin = 0.1*AliL3Transform::GetSolenoidField();
	68
	69	// Init ITS geometry
	70	fITSgeom = new AliITSgeom();
	71	fITSgeom->ReadNewFile("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymmFMD.det");
	72	if (!fITSgeom)
	73	AliError("ITS geometry not found!");
	74
c27411e7	75	// Init PID
c27411e7	76	AliPID pid;
e10815f1	77	}
	78
	79	//_____________________________________________________________________________
	80	Bool_t AliHoughFilter::Filter(AliRawEvent* event, AliESD* esd)
	81	{
d6186083	82	// Run fast online reconstruction
	83	// based on the HLT tracking algorithms
	84
	85	TStopwatch globaltimer;
	86	globaltimer.Start();
	87
d6186083	88	TStopwatch timer;
	89	timer.Start();
	90
	91	TTree *treeITSclusters = new TTree("TreeL3ITSclusters"," "); //make a tree
79eedd0f	92	treeITSclusters->SetDirectory(0);
79eedd0f	93
d6186083	94	RunITSclusterer(event,treeITSclusters);
	95	RunITSvertexer(esd,treeITSclusters);
	96	RunTPCtracking(event,esd);
	97	RunITStracking(esd,treeITSclusters);
	98	delete treeITSclusters;
	99
d6186083	100	AliInfo(Form("Event filter has finished in %f seconds\n\n\n\n\n\n",globaltimer.RealTime()));
	101
	102	return kTRUE;
	103	}
	104
	105	//_____________________________________________________________________________
	106	void AliHoughFilter::RunITSclusterer(AliRawEvent* event, TTree *treeClusters)
	107	{
	108	// Run ITS Clusterer
	109	// The clusters are stored in a tree
	110
	111	TStopwatch timer;
	112	timer.Start();
	113
	114	if(!fITSgeom)
	115	AliError("ITS geometry not created!");
	116	AliL3ITSclusterer clusterer(fITSgeom);
	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(AliESD* esd, TTree *treeClusters)
	127	{
	128	// Run SPD vertexerZ
	129	// Store the result in the ESD
	130
	131	TStopwatch timer;
	132	timer.Start();
e10815f1	133
d6186083	134	AliL3ITSVertexerZ 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(AliRawEvent* event, AliESD* 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
e10815f1	147	TStopwatch timer;
	148	timer.Start();
	149
d6186083	150	const AliESDVertex *vertex = esd->GetVertex();
d6186083	151	Float_t zvertex = vertex->GetZv();
e10815f1	152
	153	AliL3Hough *hough1 = new AliL3Hough();
	154
	155	hough1->SetThreshold(4);
d6186083	156	hough1->CalcTransformerParams(fPtmin);
e10815f1	157	hough1->SetPeakThreshold(70,-1);
d6186083	158	hough1->Init(100,4,event,zvertex);
e10815f1	159	hough1->SetAddHistograms();
	160
	161	AliL3Hough *hough2 = new AliL3Hough();
d6186083	162
e10815f1	163	hough2->SetThreshold(4);
d6186083	164	hough2->CalcTransformerParams(fPtmin);
e10815f1	165	hough2->SetPeakThreshold(70,-1);
d6186083	166	hough2->Init(100,4,event,zvertex);
e10815f1	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())
d6186083	175	::Fatal("AliL3Hough::WaitForThreadFinish"," Can not join the required thread! ");
e10815f1	176	if(hough2->WaitForThreadFinish())
d6186083	177	::Fatal("AliL3Hough::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	*/
e10815f1	197
	198	nglobaltracks += hough1->FillESD(esd);
	199	nglobaltracks += hough2->FillESD(esd);
	200
d6186083	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	*/
e10815f1	207
	208	delete hough1;
	209	delete hough2;
	210
d6186083	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(AliESD* 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	AliL3ITStracker itsTracker(fITSgeom);
	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()));
e10815f1	236
e10815f1	237	}