[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
  AliKalmanTrack::SetConvConst(
     1000/0.299792458/AliL3Transform::GetSolenoidField()
  );
  AliMagF* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
  AliTracker::SetFieldMap(field);
  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!");

}

//_____________________________________________________________________________
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
	64	// Init magnetic field
	65	AliKalmanTrack::SetConvConst(
	66	1000/0.299792458/AliL3Transform::GetSolenoidField()
	67	);
	68	AliMagF* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
	69	AliTracker::SetFieldMap(field);
	70	fPtmin = 0.1*AliL3Transform::GetSolenoidField();
	71
	72	// Init ITS geometry
	73	fITSgeom = new AliITSgeom();
	74	fITSgeom->ReadNewFile("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymmFMD.det");
	75	if (!fITSgeom)
	76	AliError("ITS geometry not found!");
	77
e10815f1	78	}
	79
	80	//_____________________________________________________________________________
	81	Bool_t AliHoughFilter::Filter(AliRawEvent* event, AliESD* esd)
	82	{
d6186083	83	// Run fast online reconstruction
	84	// based on the HLT tracking algorithms
	85
	86	TStopwatch globaltimer;
	87	globaltimer.Start();
	88
d6186083	89	TStopwatch timer;
	90	timer.Start();
	91
	92	TTree *treeITSclusters = new TTree("TreeL3ITSclusters"," "); //make a tree
79eedd0f	93	treeITSclusters->SetDirectory(0);
79eedd0f	94
d6186083	95	RunITSclusterer(event,treeITSclusters);
	96	RunITSvertexer(esd,treeITSclusters);
	97	RunTPCtracking(event,esd);
	98	RunITStracking(esd,treeITSclusters);
	99	delete treeITSclusters;
	100
d6186083	101	AliInfo(Form("Event filter has finished in %f seconds\n\n\n\n\n\n",globaltimer.RealTime()));
	102
	103	return kTRUE;
	104	}
	105
	106	//_____________________________________________________________________________
	107	void AliHoughFilter::RunITSclusterer(AliRawEvent* event, TTree *treeClusters)
	108	{
	109	// Run ITS Clusterer
	110	// The clusters are stored in a tree
	111
	112	TStopwatch timer;
	113	timer.Start();
	114
	115	if(!fITSgeom)
	116	AliError("ITS geometry not created!");
	117	AliL3ITSclusterer clusterer(fITSgeom);
	118	AliRawReader *itsrawreader=new AliRawReaderRoot(event);
	119	clusterer.Digits2Clusters(itsrawreader,treeClusters);
	120	delete itsrawreader;
	121	AliInfo(Form("ITS clusterer has finished in %f seconds\n",timer.RealTime()));
	122
	123	}
	124
	125
	126	//_____________________________________________________________________________
	127	void AliHoughFilter::RunITSvertexer(AliESD* esd, TTree *treeClusters)
	128	{
	129	// Run SPD vertexerZ
	130	// Store the result in the ESD
	131
	132	TStopwatch timer;
	133	timer.Start();
e10815f1	134
d6186083	135	AliL3ITSVertexerZ vertexer;
	136	AliESDVertex *vertex = vertexer.FindVertexForCurrentEvent(fITSgeom,treeClusters);
	137	esd->SetVertex(vertex);
	138	AliInfo(Form("ITS vertexer has finished in %f seconds\n",timer.RealTime()));
	139
	140	}
	141
	142	//_____________________________________________________________________________
	143	void AliHoughFilter::RunTPCtracking(AliRawEvent* event, AliESD* esd)
	144	{
	145	// Run hough transform tracking in TPC
	146	// The z of the vertex is taken from the ESD
	147	// The result of the tracking is stored in the ESD
e10815f1	148	TStopwatch timer;
	149	timer.Start();
	150
d6186083	151	const AliESDVertex *vertex = esd->GetVertex();
d6186083	152	Float_t zvertex = vertex->GetZv();
e10815f1	153
	154	AliL3Hough *hough1 = new AliL3Hough();
	155
	156	hough1->SetThreshold(4);
d6186083	157	hough1->CalcTransformerParams(fPtmin);
e10815f1	158	hough1->SetPeakThreshold(70,-1);
d6186083	159	hough1->Init(100,4,event,zvertex);
e10815f1	160	hough1->SetAddHistograms();
	161
	162	AliL3Hough *hough2 = new AliL3Hough();
d6186083	163
e10815f1	164	hough2->SetThreshold(4);
d6186083	165	hough2->CalcTransformerParams(fPtmin);
e10815f1	166	hough2->SetPeakThreshold(70,-1);
d6186083	167	hough2->Init(100,4,event,zvertex);
e10815f1	168	hough2->SetAddHistograms();
	169
	170	Int_t nglobaltracks = 0;
	171	/* In case we run HLT code in 2 threads */
	172	hough1->StartProcessInThread(0,17);
	173	hough2->StartProcessInThread(18,35);
	174
	175	if(hough1->WaitForThreadFinish())
d6186083	176	::Fatal("AliL3Hough::WaitForThreadFinish"," Can not join the required thread! ");
e10815f1	177	if(hough2->WaitForThreadFinish())
d6186083	178	::Fatal("AliL3Hough::WaitForThreadFinish"," Can not join the required thread! ");
	179
	180	/* In case we run HLT code in the main thread
	181	for(Int_t slice=0; slice<=17; slice++)
	182	{
	183	hough1->ReadData(slice,0);
	184	hough1->Transform();
	185	hough1->AddAllHistogramsRows();
	186	hough1->FindTrackCandidatesRow();
	187	hough1->AddTracks();
	188	}
	189	for(Int_t slice=18; slice<=35; slice++)
	190	{
	191	hough2->ReadData(slice,0);
	192	hough2->Transform();
	193	hough2->AddAllHistogramsRows();
	194	hough2->FindTrackCandidatesRow();
	195	hough2->AddTracks();
	196	}
	197	*/
e10815f1	198
	199	nglobaltracks += hough1->FillESD(esd);
	200	nglobaltracks += hough2->FillESD(esd);
	201
d6186083	202	/* In case we want to debug the ESD
	203	gSystem->MakeDirectory("hough1");
	204	hough1->WriteTracks("./hough1");
	205	gSystem->MakeDirectory("hough2");
	206	hough2->WriteTracks("./hough2");
	207	*/
e10815f1	208
	209	delete hough1;
	210	delete hough2;
	211
d6186083	212	AliInfo(Form("\nHough Transformer has found %d TPC tracks in %f seconds\n\n", nglobaltracks,timer.RealTime()));
	213
	214	}
	215
	216	//_____________________________________________________________________________
	217	void AliHoughFilter::RunITStracking(AliESD* esd, TTree *treeClusters)
	218	{
	219	// Run the ITS tracker
	220	// The tracks from the HT TPC tracking are used as seeds
	221	// The prologated tracks are updated in the ESD
	222	TStopwatch timer;
	223	timer.Start();
	224
	225	Double_t vtxPos[3];
	226	Double_t vtxErr[3]={0.005,0.005,0.010};
	227	const AliESDVertex *vertex = esd->GetVertex();
	228	vertex->GetXYZ(vtxPos);
	229
	230	AliL3ITStracker itsTracker(fITSgeom);
	231	itsTracker.SetVertex(vtxPos,vtxErr);
	232
	233	itsTracker.LoadClusters(treeClusters);
	234	itsTracker.Clusters2Tracks(esd);
	235	itsTracker.UnloadClusters();
	236	AliInfo(Form("ITS tracker has finished in %f seconds\n",timer.RealTime()));
e10815f1	237
e10815f1	238	}