Adding support for raw data events without subevent header
[u/mrichter/AliRoot.git] / RAW / AliHoughFilter.cxx
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 "AliL3StandardIncludes.h"
27 #include "AliL3Logging.h"
28 #include "AliL3Transform.h"
29 #include "AliL3Hough.h"
30 #include "AliLog.h"
31 #include <AliL3ITSclusterer.h>
32 #include <AliL3ITSVertexerZ.h>
33 #include <AliL3ITStracker.h>
34
35 #include "AliHoughFilter.h"
36
37 #include "AliRawReaderRoot.h"
38 #include <AliMagF.h>
39 #include <AliMagFMaps.h>
40 #include <AliKalmanTrack.h>
41 #include <AliITSgeom.h>
42
43 ClassImp(AliHoughFilter)
44
45 //_____________________________________________________________________________
46 AliHoughFilter::AliHoughFilter()
47 {
48 // default constructor
49
50   // Init debug level
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;
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))
62     AliError("HLT initialization failed!");
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
78 }
79
80 //_____________________________________________________________________________
81 Bool_t AliHoughFilter::Filter(AliRawEvent* event, AliESD* esd)
82 {
83   // Run fast online reconstruction
84   // based on the HLT tracking algorithms
85
86   TStopwatch globaltimer;
87   globaltimer.Start();
88
89   TStopwatch timer;
90   timer.Start();
91
92   TTree *treeITSclusters = new TTree("TreeL3ITSclusters"," "); //make a tree
93   treeITSclusters->SetDirectory(0);
94
95   RunITSclusterer(event,treeITSclusters);
96   RunITSvertexer(esd,treeITSclusters);
97   RunTPCtracking(event,esd);
98   RunITStracking(esd,treeITSclusters);
99   delete treeITSclusters;
100
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();
134
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
148   TStopwatch timer;
149   timer.Start();
150
151   const AliESDVertex *vertex = esd->GetVertex();
152   Float_t zvertex = vertex->GetZv();
153
154   AliL3Hough *hough1 = new AliL3Hough();
155     
156   hough1->SetThreshold(4);
157   hough1->CalcTransformerParams(fPtmin);
158   hough1->SetPeakThreshold(70,-1);
159   hough1->Init(100,4,event,zvertex);
160   hough1->SetAddHistograms();
161
162   AliL3Hough *hough2 = new AliL3Hough();
163   
164   hough2->SetThreshold(4);
165   hough2->CalcTransformerParams(fPtmin);
166   hough2->SetPeakThreshold(70,-1);
167   hough2->Init(100,4,event,zvertex);
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())
176     ::Fatal("AliL3Hough::WaitForThreadFinish"," Can not join the required thread! ");
177   if(hough2->WaitForThreadFinish())
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     */
198
199   nglobaltracks += hough1->FillESD(esd);
200   nglobaltracks += hough2->FillESD(esd);
201
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     */
208
209   delete hough1;
210   delete hough2;
211
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()));
237
238 }