[u/mrichter/AliRoot.git] / HLT / src / AliHLTReconstructor.cxx

// $Id$

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// class for HLT reconstruction                                              //
// <Cvetan.Cheshkov@cern.ch>                                                 //
///////////////////////////////////////////////////////////////////////////////

// very ugly but it has to work fast
#ifdef use_reconstruction

#include <Riostream.h>
#include <TSystem.h>
#include <TArrayF.h>

#include "AliL3StandardIncludes.h"
#include "AliL3Logging.h"
#include "AliLevel3.h"
#include "AliL3Evaluate.h"
#include "AliHLTReconstructor.h"
#include "AliL3Transform.h"
#include "AliL3Hough.h"
#include "AliL3FileHandler.h"
#include "AliL3Track.h"
#include "AliL3HoughTrack.h"
#include "AliL3TrackArray.h"
#include "AliRunLoader.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliESD.h"
#include "AliESDHLTtrack.h"

#if __GNUC__== 3
using namespace std;
#endif


ClassImp(AliHLTReconstructor)

void AliHLTReconstructor::Reconstruct(AliRunLoader* runLoader) const
{
  if(!runLoader) {
    LOG(AliL3Log::kFatal,"AliHLTReconstructor::Reconstruct","RunLoader")
      <<" Missing RunLoader! 0x0"<<ENDLOG;
  }
  gSystem->Exec("rm -rf hlt");
  gSystem->MakeDirectory("hlt");
  gSystem->Exec("rm -rf hough");
  gSystem->MakeDirectory("hough");
  gSystem->Exec("rm -f alirunfile.root");
  gSystem->Exec("ln -s galice.root alirunfile.root");

  Bool_t isinit=AliL3Transform::Init("./",kTRUE);
  if(!isinit){
    cerr << "Could not create transform settings, please check log for error messages!" << endl;
    return;
  }

  Int_t nEvents = runLoader->GetNumberOfEvents();

  for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
    runLoader->GetEvent(iEvent);

    ReconstructWithConformalMapping(runLoader,iEvent);
    ReconstructWithHoughTransform(runLoader,iEvent);
  }
}

void AliHLTReconstructor::ReconstructWithConformalMapping(AliRunLoader* runLoader,Int_t iEvent) const
{
  AliLevel3 *fHLT = new AliLevel3(runLoader);
  fHLT->Init("./", AliLevel3::kRunLoader, 1);

  Int_t phiSegments = 50;
  Int_t etaSegments = 100;
  Int_t trackletlength = 3;
  Int_t tracklength = 10;
  Int_t rowscopetracklet = 2;
  Int_t rowscopetrack = 10;
  Double_t minPtFit = 0;
  Double_t maxangle = 0.1745;
  Double_t goodDist = 5;
  Double_t maxphi = 0.1;
  Double_t maxeta = 0.1;
  Double_t hitChi2Cut = 20;
  Double_t goodHitChi2 = 5;
  Double_t trackChi2Cut = 10;
  Double_t xyerror = -1;
  Double_t zerror =  -1;
  
  fHLT->SetClusterFinderParam(xyerror,zerror,kTRUE);
  fHLT->SetTrackerParam(phiSegments, etaSegments, 
			trackletlength, tracklength,
			rowscopetracklet, rowscopetrack,
			minPtFit, maxangle, goodDist, hitChi2Cut,
			goodHitChi2, trackChi2Cut, 50, maxphi, maxeta, kTRUE);
  fHLT->SetMergerParameters(2,3,0.003,0.1,0.05);
  fHLT->DoMc();
  fHLT->WriteFiles("./hlt/");  

  fHLT->ProcessEvent(0, 35, iEvent);

  char filename[256];
  sprintf(filename, "confmap_%d",iEvent);
  fHLT->DoBench(filename);

  delete fHLT;

}

void AliHLTReconstructor::ReconstructWithHoughTransform(AliRunLoader* runLoader,Int_t iEvent) const
{
  Float_t ptmin = 0.1*AliL3Transform::GetSolenoidField();

  Float_t zvertex = 0;
  TArrayF mcVertex(3); 
  runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
  zvertex = mcVertex[2];

  cout<<" Hough Tranform will run with ptmin="<<ptmin<<" and zvertex="<<zvertex<<endl;

  AliL3Hough *hough = new AliL3Hough();
    
  hough->SetThreshold(4);
  hough->SetTransformerParams(76,140,ptmin,-1);
  hough->SetPeakThreshold(70,-1);
  hough->SetRunLoader(runLoader);
  hough->Init("./", kFALSE, 100, kFALSE,4,0,0,zvertex);
  hough->SetAddHistograms();

  for(int slice=0; slice<=35; slice++)
    {
      //     cout<<"Processing slice "<<slice<<endl;
      hough->ReadData(slice,iEvent);
      hough->Transform();
      hough->AddAllHistogramsRows();
      hough->FindTrackCandidatesRow();
      //     hough->WriteTracks(slice,"./hough");
      hough->AddTracks();
    }
  hough->WriteTracks("./hough");
  
  char filename[256];
  sprintf(filename, "hough_%d",iEvent);
  hough->DoBench(filename);

  delete hough;
}

void AliHLTReconstructor::FillESD(AliRunLoader* runLoader, 
				  AliESD* esd) const
{
  Int_t iEvent = runLoader->GetEventNumber();

  FillESDforConformalMapping(esd,iEvent);
  FillESDforHoughTransform(esd,iEvent);
}

void AliHLTReconstructor::FillESDforConformalMapping(AliESD* esd,Int_t iEvent) const
{
  //Assign MC labels for found tracks
  int slicerange[2]={0,35};
  int good = (int)(0.4*AliL3Transform::GetNRows());
  int nclusters = (int)(0.4*AliL3Transform::GetNRows());
  float ptmin = 0.;
  float ptmax = 0.;
  float maxfalseratio = 0.1;
  
  AliL3Evaluate *fHLTEval = new AliL3Evaluate("./hlt",nclusters,good,ptmin,ptmax,slicerange);
  fHLTEval->SetMaxFalseClusters(maxfalseratio);

  fHLTEval->LoadData(iEvent,kTRUE);
  fHLTEval->AssignPIDs();
  fHLTEval->AssignIDs();
  AliL3TrackArray *fTracks = fHLTEval->GetTracks();
  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliL3Track *tpt = (AliL3Track *)fTracks->GetCheckedTrack(i);
      if(!tpt) continue; 
      
      AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ; 

      esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
      esdtrack->SetNHits(tpt->GetNHits());
      esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
      esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
      esdtrack->SetPt(tpt->GetPt());
      esdtrack->SetPsi(tpt->GetPsi());
      esdtrack->SetTgl(tpt->GetTgl());
      esdtrack->SetCharge(tpt->GetCharge());
      esdtrack->SetMCid(tpt->GetMCid());
      esdtrack->SetSector(tpt->GetSector());
      esdtrack->SetPID(tpt->GetPID());
      esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());

      esd->AddHLTConfMapTrack(esdtrack);
      delete esdtrack;
    }

  delete fHLTEval;
}

void AliHLTReconstructor::FillESDforHoughTransform(AliESD* esd,Int_t iEvent) const
{
  char filename[256];
  sprintf(filename,"./hough/tracks_%d.raw",iEvent);
  
  AliL3FileHandler *tfile = new AliL3FileHandler();
  if(!tfile->SetBinaryInput(filename)){
    Error("FillESD","Inputfile ",filename," does not exist");
    return;
  }
  
  AliL3TrackArray *fTracks = new AliL3TrackArray("AliL3HoughTrack");
  tfile->Binary2TrackArray(fTracks);
  tfile->CloseBinaryInput();
  delete tfile;
  
  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliL3HoughTrack *tpt = (AliL3HoughTrack *)fTracks->GetCheckedTrack(i);
      if(!tpt) continue; 
      
      AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ; 

      esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
      esdtrack->SetNHits(tpt->GetNHits());
      esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
      esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
      esdtrack->SetPt(tpt->GetPt());
      esdtrack->SetPsi(tpt->GetPsi());
      esdtrack->SetTgl(tpt->GetTgl());
      esdtrack->SetCharge(tpt->GetCharge());
      esdtrack->SetMCid(tpt->GetMCid());
      esdtrack->SetWeight(tpt->GetWeight());
      esdtrack->SetSector(tpt->GetSector());
      esdtrack->SetBinXY(tpt->GetBinX(),tpt->GetBinY(),tpt->GetSizeX(),tpt->GetSizeY());
      esdtrack->SetPID(tpt->GetPID());
      esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());

      esd->AddHLTHoughTrack(esdtrack);
      delete esdtrack;
    }

  delete fTracks;
}

#endif
Commit	Line	Data
de3c3890	1	// $Id$
	2
	3	///////////////////////////////////////////////////////////////////////////////
	4	// //
	5	// class for HLT reconstruction //
	6	// <Cvetan.Cheshkov@cern.ch> //
	7	///////////////////////////////////////////////////////////////////////////////
	8
	9	// very ugly but it has to work fast
	10	#ifdef use_reconstruction
	11
	12	#include <Riostream.h>
	13	#include <TSystem.h>
	14	#include <TArrayF.h>
	15
	16	#include "AliL3StandardIncludes.h"
9a74ed1c	17	#include "AliL3Logging.h"
de3c3890	18	#include "AliLevel3.h"
	19	#include "AliL3Evaluate.h"
	20	#include "AliHLTReconstructor.h"
	21	#include "AliL3Transform.h"
	22	#include "AliL3Hough.h"
	23	#include "AliL3FileHandler.h"
	24	#include "AliL3Track.h"
	25	#include "AliL3HoughTrack.h"
	26	#include "AliL3TrackArray.h"
	27	#include "AliRunLoader.h"
	28	#include "AliHeader.h"
	29	#include "AliGenEventHeader.h"
	30	#include "AliESD.h"
	31	#include "AliESDHLTtrack.h"
	32
	33	#if __GNUC__== 3
	34	using namespace std;
	35	#endif
	36
	37
	38	ClassImp(AliHLTReconstructor)
	39
	40	void AliHLTReconstructor::Reconstruct(AliRunLoader* runLoader) const
	41	{
9a74ed1c	42	if(!runLoader) {
	43	LOG(AliL3Log::kFatal,"AliHLTReconstructor::Reconstruct","RunLoader")
	44	<<" Missing RunLoader! 0x0"<<ENDLOG;
	45	}
de3c3890	46	gSystem->Exec("rm -rf hlt");
	47	gSystem->MakeDirectory("hlt");
	48	gSystem->Exec("rm -rf hough");
	49	gSystem->MakeDirectory("hough");
	50	gSystem->Exec("rm -f alirunfile.root");
	51	gSystem->Exec("ln -s galice.root alirunfile.root");
	52
	53	Bool_t isinit=AliL3Transform::Init("./",kTRUE);
	54	if(!isinit){
	55	cerr << "Could not create transform settings, please check log for error messages!" << endl;
	56	return;
	57	}
	58
	59	Int_t nEvents = runLoader->GetNumberOfEvents();
	60
	61	for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
	62	runLoader->GetEvent(iEvent);
	63
	64	ReconstructWithConformalMapping(runLoader,iEvent);
	65	ReconstructWithHoughTransform(runLoader,iEvent);
	66	}
	67	}
	68
	69	void AliHLTReconstructor::ReconstructWithConformalMapping(AliRunLoader* runLoader,Int_t iEvent) const
	70	{
	71	AliLevel3 *fHLT = new AliLevel3(runLoader);
	72	fHLT->Init("./", AliLevel3::kRunLoader, 1);
	73
	74	Int_t phiSegments = 50;
	75	Int_t etaSegments = 100;
	76	Int_t trackletlength = 3;
	77	Int_t tracklength = 10;
	78	Int_t rowscopetracklet = 2;
	79	Int_t rowscopetrack = 10;
	80	Double_t minPtFit = 0;
	81	Double_t maxangle = 0.1745;
	82	Double_t goodDist = 5;
	83	Double_t maxphi = 0.1;
	84	Double_t maxeta = 0.1;
	85	Double_t hitChi2Cut = 20;
	86	Double_t goodHitChi2 = 5;
	87	Double_t trackChi2Cut = 10;
	88	Double_t xyerror = -1;
	89	Double_t zerror = -1;
	90
	91	fHLT->SetClusterFinderParam(xyerror,zerror,kTRUE);
	92	fHLT->SetTrackerParam(phiSegments, etaSegments,
	93	trackletlength, tracklength,
	94	rowscopetracklet, rowscopetrack,
	95	minPtFit, maxangle, goodDist, hitChi2Cut,
	96	goodHitChi2, trackChi2Cut, 50, maxphi, maxeta, kTRUE);
	97	fHLT->SetMergerParameters(2,3,0.003,0.1,0.05);
	98	fHLT->DoMc();
	99	fHLT->WriteFiles("./hlt/");
	100
	101	fHLT->ProcessEvent(0, 35, iEvent);
	102
	103	char filename[256];
	104	sprintf(filename, "confmap_%d",iEvent);
	105	fHLT->DoBench(filename);
	106
	107	delete fHLT;
	108
	109	}
110
111	void AliHLTReconstructor::ReconstructWithHoughTransform(AliRunLoader* runLoader,Int_t iEvent) const
112	{
113	Float_t ptmin = 0.1*AliL3Transform::GetSolenoidField();
114
115	Float_t zvertex = 0;
116	TArrayF mcVertex(3);
117	runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
118	zvertex = mcVertex[2];
119
120	cout<<" Hough Tranform will run with ptmin="<<ptmin<<" and zvertex="<<zvertex<<endl;
121
122	AliL3Hough *hough = new AliL3Hough();
123
124	hough->SetThreshold(4);
125	hough->SetTransformerParams(76,140,ptmin,-1);
126	hough->SetPeakThreshold(70,-1);
127	hough->SetRunLoader(runLoader);
128	hough->Init("./", kFALSE, 100, kFALSE,4,0,0,zvertex);
129	hough->SetAddHistograms();
130
131	for(int slice=0; slice<=35; slice++)
132	{
133	// cout<<"Processing slice "<<slice<<endl;
134	hough->ReadData(slice,iEvent);
135	hough->Transform();
136	hough->AddAllHistogramsRows();
137	hough->FindTrackCandidatesRow();
138	// hough->WriteTracks(slice,"./hough");
139	hough->AddTracks();
140	}
141	hough->WriteTracks("./hough");
142
143	char filename[256];
144	sprintf(filename, "hough_%d",iEvent);
145	hough->DoBench(filename);
146
147	delete hough;
148	}
149
150	void AliHLTReconstructor::FillESD(AliRunLoader* runLoader,
151	AliESD* esd) const
152	{
153	Int_t iEvent = runLoader->GetEventNumber();
154
155	FillESDforConformalMapping(esd,iEvent);
156	FillESDforHoughTransform(esd,iEvent);
157	}
158
159	void AliHLTReconstructor::FillESDforConformalMapping(AliESD* esd,Int_t iEvent) const
160	{
161	//Assign MC labels for found tracks
162	int slicerange[2]={0,35};
163	int good = (int)(0.4*AliL3Transform::GetNRows());
164	int nclusters = (int)(0.4*AliL3Transform::GetNRows());
165	float ptmin = 0.;
166	float ptmax = 0.;
167	float maxfalseratio = 0.1;
168
169	AliL3Evaluate *fHLTEval = new AliL3Evaluate("./hlt",nclusters,good,ptmin,ptmax,slicerange);
170	fHLTEval->SetMaxFalseClusters(maxfalseratio);
171
172	fHLTEval->LoadData(iEvent,kTRUE);
173	fHLTEval->AssignPIDs();
174	fHLTEval->AssignIDs();
175	AliL3TrackArray *fTracks = fHLTEval->GetTracks();
176	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
177	{
178	AliL3Track tpt = (AliL3Track )fTracks->GetCheckedTrack(i);
179	if(!tpt) continue;
180
181	AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ;
182
183	esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
184	esdtrack->SetNHits(tpt->GetNHits());
185	esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
186	esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
187	esdtrack->SetPt(tpt->GetPt());
188	esdtrack->SetPsi(tpt->GetPsi());
189	esdtrack->SetTgl(tpt->GetTgl());
190	esdtrack->SetCharge(tpt->GetCharge());
191	esdtrack->SetMCid(tpt->GetMCid());
192	esdtrack->SetSector(tpt->GetSector());
193	esdtrack->SetPID(tpt->GetPID());
194	esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());
195
196	esd->AddHLTConfMapTrack(esdtrack);
197	delete esdtrack;
198	}
199
200	delete fHLTEval;
201	}
202
203	void AliHLTReconstructor::FillESDforHoughTransform(AliESD* esd,Int_t iEvent) const
204	{
205	char filename[256];
206	sprintf(filename,"./hough/tracks_%d.raw",iEvent);
207
208	AliL3FileHandler *tfile = new AliL3FileHandler();
209	if(!tfile->SetBinaryInput(filename)){
210	Error("FillESD","Inputfile ",filename," does not exist");
211	return;
212	}
213
214	AliL3TrackArray *fTracks = new AliL3TrackArray("AliL3HoughTrack");
215	tfile->Binary2TrackArray(fTracks);
216	tfile->CloseBinaryInput();
217	delete tfile;
218
219	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
220	{
221	AliL3HoughTrack tpt = (AliL3HoughTrack )fTracks->GetCheckedTrack(i);
222	if(!tpt) continue;
223
224	AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ;
225
226	esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
227	esdtrack->SetNHits(tpt->GetNHits());
228	esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
229	esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
230	esdtrack->SetPt(tpt->GetPt());
231	esdtrack->SetPsi(tpt->GetPsi());
232	esdtrack->SetTgl(tpt->GetTgl());
233	esdtrack->SetCharge(tpt->GetCharge());
234	esdtrack->SetMCid(tpt->GetMCid());
235	esdtrack->SetWeight(tpt->GetWeight());
236	esdtrack->SetSector(tpt->GetSector());
237	esdtrack->SetBinXY(tpt->GetBinX(),tpt->GetBinY(),tpt->GetSizeX(),tpt->GetSizeY());
238	esdtrack->SetPID(tpt->GetPID());
239	esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());
240
241	esd->AddHLTHoughTrack(esdtrack);
242	delete esdtrack;
243	}
244
245	delete fTracks;
246	}
247
248	#endif