Changes to integrate the MLUC libs from
[u/mrichter/AliRoot.git] / HLT / src / AliHLTReconstructor.cxx
1 // $Id$
2
3 ///////////////////////////////////////////////////////////////////////////////
4 //                                                                           //
5 // class for HLT reconstruction                                              //
6 // <Cvetan.Cheshkov@cern.ch>                                                 //
7 // <loizides@ikf.uni-frankfurt.de>                                           //
8 ///////////////////////////////////////////////////////////////////////////////
9
10 // very ugly but it has to work fast
11 #ifdef use_reconstruction
12
13 #include <Riostream.h>
14 #include <TSystem.h>
15 #include <TArrayF.h>
16
17 #include <AliRunLoader.h>
18 #include <AliHeader.h>
19 #include <AliGenEventHeader.h>
20 #include <AliESD.h>
21 #include <AliESDHLTtrack.h>
22
23 #include "AliL3StandardIncludes.h"
24 #include "AliL3Logging.h"
25 #include "AliLevel3.h"
26 #include "AliL3Evaluate.h"
27 #include "AliHLTReconstructor.h"
28 #include "AliL3Transform.h"
29 #include "AliL3Hough.h"
30 #include "AliL3FileHandler.h"
31 #include "AliL3Track.h"
32 #include "AliL3HoughTrack.h"
33 #include "AliL3TrackArray.h"
34
35 #if __GNUC__== 3
36 using namespace std;
37 #endif
38
39
40 ClassImp(AliHLTReconstructor)
41
42 AliHLTReconstructor::AliHLTReconstructor(): AliReconstructor() 
43
44   //constructor
45 #ifndef use_logging
46   AliL3Log::fgLevel=AliL3Log::kWarning;
47 #endif
48   fDoTracker=1;
49   fDoHough=1;
50   fDoBench=0;
51   fDoCleanUp=1;
52 }
53
54 AliHLTReconstructor::AliHLTReconstructor(Bool_t doTracker, Bool_t doHough): AliReconstructor() 
55
56   //constructor
57 #ifndef use_logging
58   AliL3Log::fgLevel=AliL3Log::kWarning;
59 #endif
60   fDoTracker=doTracker;
61   fDoHough=doHough;
62   fDoBench=0;
63   fDoCleanUp=1;
64 }
65
66 AliHLTReconstructor::~AliHLTReconstructor()
67
68   //deconstructor
69   if(fDoCleanUp){
70     char name[256];
71     gSystem->Exec("rm -rf hlt");
72     sprintf(name, "rm -f confmap_*.root confmap_*.dat");
73     gSystem->Exec(name);
74     gSystem->Exec("rm -rf hough");
75     sprintf(name, "rm -f hough_*.root hough_*.dat");
76     gSystem->Exec(name);
77   }
78 }
79
80 void AliHLTReconstructor::Reconstruct(AliRunLoader* runLoader) const
81 {
82   // do the standard and hough reconstruction chain
83   if(!runLoader) {
84     LOG(AliL3Log::kFatal,"AliHLTReconstructor::Reconstruct","RunLoader")
85       <<" Missing RunLoader! 0x0"<<ENDLOG;
86     return;
87   }
88   gSystem->Exec("rm -rf hlt");
89   gSystem->MakeDirectory("hlt");
90   gSystem->Exec("rm -rf hough");
91   gSystem->MakeDirectory("hough");
92
93   Bool_t isinit=AliL3Transform::Init(runLoader);
94   if(!isinit){
95     LOG(AliL3Log::kError,"AliHLTReconstructor::Reconstruct","Transformer")
96      << "Could not create transform settings, please check log for error messages!" << ENDLOG;
97     return;
98   }
99
100   Int_t nEvents = runLoader->GetNumberOfEvents();
101
102   for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
103     runLoader->GetEvent(iEvent);
104
105     if(fDoTracker) ReconstructWithConformalMapping(runLoader,iEvent);
106     if(fDoHough) ReconstructWithHoughTransform(runLoader,iEvent);
107   }
108 }
109
110 void AliHLTReconstructor::ReconstructWithConformalMapping(AliRunLoader* runLoader,Int_t iEvent) const
111 {
112   // reconstruct with conformal mapper
113   AliLevel3 *fHLT = new AliLevel3(runLoader);
114   fHLT->Init("./", AliLevel3::kRunLoader, 1);
115
116   Int_t phiSegments = 50;
117   Int_t etaSegments = 100;
118   Int_t trackletlength = 3;
119   Int_t tracklength = 10;
120   Int_t rowscopetracklet = 2;
121   Int_t rowscopetrack = 10;
122   Double_t minPtFit = 0;
123   Double_t maxangle = 0.1745;
124   Double_t goodDist = 5;
125   Double_t maxphi = 0.1;
126   Double_t maxeta = 0.1;
127   Double_t hitChi2Cut = 20;
128   Double_t goodHitChi2 = 5;
129   Double_t trackChi2Cut = 10;
130   Double_t xyerror = -1;
131   Double_t zerror =  -1;
132   
133   fHLT->SetClusterFinderParam(xyerror,zerror,kTRUE);
134   fHLT->SetTrackerParam(phiSegments, etaSegments, 
135                         trackletlength, tracklength,
136                         rowscopetracklet, rowscopetrack,
137                         minPtFit, maxangle, goodDist, hitChi2Cut,
138                         goodHitChi2, trackChi2Cut, 50, maxphi, maxeta, kTRUE);
139   fHLT->SetTrackerParam(phiSegments, etaSegments, 
140                         trackletlength, tracklength,
141                         rowscopetracklet, rowscopetrack,
142                         minPtFit, maxangle, goodDist, hitChi2Cut,
143                         goodHitChi2, trackChi2Cut, 50, maxphi, maxeta, kFALSE);
144   fHLT->SetMergerParameters(2,3,0.003,0.1,0.05);
145   fHLT->DoMc();
146   fHLT->DoNonVertexTracking(); /*2 tracking passes, last without vertex contraint.*/
147   fHLT->WriteFiles("./hlt/");  
148   fHLT->ProcessEvent(0, 35, iEvent);
149   if(fDoBench){
150     char filename[256];
151     sprintf(filename, "confmap_%d",iEvent);
152     fHLT->DoBench(filename);
153   }
154
155   delete fHLT;
156 }
157
158 void AliHLTReconstructor::ReconstructWithHoughTransform(AliRunLoader* runLoader,Int_t iEvent) const
159 {
160   //reconstruct with hough
161   Float_t ptmin = 0.1*AliL3Transform::GetSolenoidField();
162
163   Float_t zvertex = 0;
164   TArrayF mcVertex(3); 
165   AliHeader * header = runLoader->GetHeader();
166   if (header) {
167     AliGenEventHeader * genHeader = header->GenEventHeader();
168     if (genHeader) genHeader->PrimaryVertex(mcVertex);
169   }
170   zvertex = mcVertex[2];
171
172   LOG(AliL3Log::kInformational,"AliHLTReconstructor::Reconstruct","HoughTransform")
173     <<" Hough Transform will run with ptmin="<<ptmin<<" and zvertex="<<zvertex<<ENDLOG;
174
175   AliL3Hough *hough = new AliL3Hough();
176     
177   hough->SetThreshold(4);
178   hough->SetTransformerParams(76,140,ptmin,-1);
179   hough->SetPeakThreshold(70,-1);
180   hough->SetRunLoader(runLoader);
181   hough->Init("./", kFALSE, 100, kFALSE,4,0,0,zvertex);
182   hough->SetAddHistograms();
183
184   for(Int_t slice=0; slice<=35; slice++)
185     {
186       //cout<<"Processing slice "<<slice<<endl;
187       hough->ReadData(slice,iEvent);
188       hough->Transform();
189       hough->AddAllHistogramsRows();
190       hough->FindTrackCandidatesRow();
191       //hough->WriteTracks(slice,"./hough");
192       hough->AddTracks();
193     }
194   hough->WriteTracks("./hough");
195   
196   if(fDoBench){
197     char filename[256];
198     sprintf(filename, "hough_%d",iEvent);
199     hough->DoBench(filename);
200   }
201   delete hough;
202 }
203
204 void AliHLTReconstructor::FillESD(AliRunLoader* runLoader, 
205                                   AliESD* esd) const
206 {
207   //fill the esd file with found tracks
208   Int_t iEvent = runLoader->GetEventNumber();
209
210   if(fDoTracker) FillESDforConformalMapping(esd,iEvent);
211   if(fDoHough) FillESDforHoughTransform(esd,iEvent);
212 }
213
214 void AliHLTReconstructor::FillESDforConformalMapping(AliESD* esd,Int_t iEvent) const
215 {
216   //fill esd with tracks from conformal mapping
217   Int_t slicerange[2]={0,35};
218   Int_t good = (int)(0.4*AliL3Transform::GetNRows());
219   Int_t nclusters = (int)(0.4*AliL3Transform::GetNRows());
220   Int_t nminpointsontracks = (int)(0.3*AliL3Transform::GetNRows());
221   Float_t ptmin = 0.;
222   Float_t ptmax = 0.;
223   Float_t maxfalseratio = 0.1;
224   
225   AliL3Evaluate *fHLTEval = new AliL3Evaluate("./hlt",nclusters,good,ptmin,ptmax,slicerange);
226   fHLTEval->SetMaxFalseClusters(maxfalseratio);
227   fHLTEval->LoadData(iEvent,kTRUE);
228   fHLTEval->AssignPIDs();
229   fHLTEval->AssignIDs();
230   AliL3TrackArray *fTracks = fHLTEval->GetTracks();
231   if(!fTracks){
232     delete fHLTEval;
233     return;
234   }
235   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
236     {
237       AliL3Track *tpt = (AliL3Track *)fTracks->GetCheckedTrack(i);
238       if(!tpt) continue; 
239       if(tpt->GetNumberOfPoints() < nminpointsontracks) continue;
240       
241       AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ; 
242
243       esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
244       esdtrack->SetNHits(tpt->GetNHits());
245       esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
246       esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
247       esdtrack->SetPt(tpt->GetPt());
248       esdtrack->SetPsi(tpt->GetPsi());
249       esdtrack->SetTgl(tpt->GetTgl());
250       esdtrack->SetCharge(tpt->GetCharge());
251       esdtrack->SetMCid(tpt->GetMCid());
252       esdtrack->SetSector(tpt->GetSector());
253       esdtrack->SetPID(tpt->GetPID());
254       esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());
255
256       esd->AddHLTConfMapTrack(esdtrack);
257       delete esdtrack;
258     }
259   delete fHLTEval;
260 }
261
262 void AliHLTReconstructor::FillESDforHoughTransform(AliESD* esd,Int_t iEvent) const
263 {
264   //fill esd with tracks from hough
265   char filename[256];
266   sprintf(filename,"./hough/tracks_%d.raw",iEvent);
267   
268   AliL3FileHandler *tfile = new AliL3FileHandler();
269   if(!tfile->SetBinaryInput(filename)){
270     LOG(AliL3Log::kError,"AliHLTReconstructor::FillESDforHoughTransform","Input file")
271       <<" Missing file "<<filename<<ENDLOG;
272     return;
273   }
274   
275   AliL3TrackArray *fTracks = new AliL3TrackArray("AliL3HoughTrack");
276   tfile->Binary2TrackArray(fTracks);
277   tfile->CloseBinaryInput();
278   delete tfile;
279   if(!fTracks) return; 
280   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
281     {
282       AliL3HoughTrack *tpt = (AliL3HoughTrack *)fTracks->GetCheckedTrack(i);
283       if(!tpt) continue; 
284       
285       AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ; 
286
287       esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
288       esdtrack->SetNHits(tpt->GetNHits());
289       esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
290       esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
291       esdtrack->SetPt(tpt->GetPt());
292       esdtrack->SetPsi(tpt->GetPsi());
293       esdtrack->SetTgl(tpt->GetTgl());
294       esdtrack->SetCharge(tpt->GetCharge());
295       esdtrack->SetMCid(tpt->GetMCid());
296       esdtrack->SetWeight(tpt->GetWeight());
297       esdtrack->SetSector(tpt->GetSector());
298       esdtrack->SetBinXY(tpt->GetBinX(),tpt->GetBinY(),tpt->GetSizeX(),tpt->GetSizeY());
299       esdtrack->SetPID(tpt->GetPID());
300       esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());
301
302       esd->AddHLTHoughTrack(esdtrack);
303       delete esdtrack;
304     }
305
306   delete fTracks;
307 }
308 #endif