Correction for Anti-Merging cut - if Track Points are not available such a pair is...
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
CommitLineData
596a855f 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// class for running the reconstruction //
21// //
22// Clusters and tracks are created for all detectors and all events by //
23// typing: //
24// //
25// AliReconstruction rec; //
26// rec.Run(); //
27// //
28// The Run method returns kTRUE in case of successful execution. //
29// The name of the galice file can be changed from the default //
30// "galice.root" by //
31// //
32// rec.SetGAliceFile("..."); //
33// //
34// The reconstruction can be switched on or off for individual detectors by //
35// //
36// rec.SetRunReconstruction("..."); //
37// //
38// The argument is a (case sensitive) string with the names of the //
39// detectors separated by a space. The special string "ALL" selects all //
40// available detectors. This is the default. //
41// //
42// The tracking in ITS, TPC and TRD and the creation of ESD tracks can be //
43// switched off by //
44// //
45// rec.SetRunTracking(kFALSE); //
46// //
47// The filling of additional ESD information can be steered by //
48// //
49// rec.SetFillESD("..."); //
50// //
51// Again, the string specifies the list of detectors. The default is "ALL". //
52// //
53// The reconstruction requires digits as input. For the creation of digits //
54// have a look at the class AliSimulation. //
55// //
56///////////////////////////////////////////////////////////////////////////////
57
58
59#include "AliReconstruction.h"
60#include "AliRunLoader.h"
61#include "AliRun.h"
62#include "AliModule.h"
63#include "AliDetector.h"
64#include "AliTracker.h"
65#include "AliESD.h"
66#include "AliHeader.h"
67#include "AliGenEventHeader.h"
68#include "AliESDpid.h"
69#include <TArrayF.h>
70
71
72ClassImp(AliReconstruction)
73
74
75//_____________________________________________________________________________
76AliReconstruction::AliReconstruction(const char* name, const char* title) :
77 TNamed(name, title)
78{
79// create reconstruction object with default parameters
80
81 Init();
82}
83
84//_____________________________________________________________________________
85AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
86 TNamed(rec)
87{
88// copy constructor
89
90 fRunReconstruction = rec.fRunReconstruction;
91 fRunTracking = rec.fRunTracking;
92 fStopOnError = rec.fStopOnError;
93
94 fGAliceFileName = rec.fGAliceFileName;
95
96 fRunLoader = NULL;
97}
98
99//_____________________________________________________________________________
100AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
101{
102// assignment operator
103
104 this->~AliReconstruction();
105 new(this) AliReconstruction(rec);
106 return *this;
107}
108
109//_____________________________________________________________________________
110AliReconstruction::~AliReconstruction()
111{
112// clean up
113
114}
115
116//_____________________________________________________________________________
117void AliReconstruction::Init()
118{
119// set default parameters
120
121 fRunReconstruction = "ALL";
122 fRunTracking = kTRUE;
123 fFillESD = "ALL";
124 fStopOnError = kFALSE;
125
126 fGAliceFileName = "galice.root";
127
128 fRunLoader = NULL;
129}
130
131
132//_____________________________________________________________________________
133void AliReconstruction::SetGAliceFile(const char* fileName)
134{
135// set the name of the galice file
136
137 fGAliceFileName = fileName;
138}
139
140
141//_____________________________________________________________________________
142Bool_t AliReconstruction::Run()
143{
144// run the reconstruction
145
146 // open the run loader
147 if (fRunLoader) delete fRunLoader;
148 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
149 if (!fRunLoader) {
150 Error("Run", "no run loader found in file %s",
151 fGAliceFileName.Data());
152 return kFALSE;
153 }
154 fRunLoader->LoadgAlice();
155 gAlice = fRunLoader->GetAliRun();
156 if (!gAlice) {
157 Error("Run", "no gAlice object found in file %s",
158 fGAliceFileName.Data());
159 return kFALSE;
160 }
161
162 // local reconstruction
163 if (!fRunReconstruction.IsNull()) {
164 if (!RunReconstruction(fRunReconstruction)) {
165 if (fStopOnError) return kFALSE;
166 }
167 }
168 if (!fRunTracking && fFillESD.IsNull()) return kTRUE;
169
170 // get loaders and trackers
171 fITSLoader = fRunLoader->GetLoader("ITSLoader");
172 if (!fITSLoader) {
173 Error("Run", "no ITS loader found");
174 if (fStopOnError) return kFALSE;
175 }
176 fITSTracker = NULL;
177 if (gAlice->GetDetector("ITS")) {
178 fITSTracker = gAlice->GetDetector("ITS")->CreateTracker();
179 }
180 if (!fITSTracker) {
181 Error("Run", "couldn't create a tracker for ITS");
182 if (fStopOnError) return kFALSE;
183 }
184
185 fTPCLoader = fRunLoader->GetLoader("TPCLoader");
186 if (!fTPCLoader) {
187 Error("Run", "no TPC loader found");
188 if (fStopOnError) return kFALSE;
189 }
190 fTPCTracker = NULL;
191 if (gAlice->GetDetector("TPC")) {
192 fTPCTracker = gAlice->GetDetector("TPC")->CreateTracker();
193 }
194 if (!fTPCTracker) {
195 Error("Run", "couldn't create a tracker for TPC");
196 if (fStopOnError) return kFALSE;
197 }
198
199 fTRDLoader = fRunLoader->GetLoader("TRDLoader");
200 if (!fTRDLoader) {
201 Error("Run", "no TRD loader found");
202 if (fStopOnError) return kFALSE;
203 }
204 fTRDTracker = NULL;
205 if (gAlice->GetDetector("TRD")) {
206 fTRDTracker = gAlice->GetDetector("TRD")->CreateTracker();
207 }
208 if (!fTRDTracker) {
209 Error("Run", "couldn't create a tracker for TRD");
210 if (fStopOnError) return kFALSE;
211 }
212
213 fTOFLoader = fRunLoader->GetLoader("TOFLoader");
214 if (!fTOFLoader) {
215 Error("Run", "no TOF loader found");
216 if (fStopOnError) return kFALSE;
217 }
218 fTOFTracker = NULL;
219 if (gAlice->GetDetector("TOF")) {
220 fTOFTracker = gAlice->GetDetector("TOF")->CreateTracker();
221 }
222 if (!fTOFTracker) {
223 Error("Run", "couldn't create a tracker for TOF");
224 if (fStopOnError) return kFALSE;
225 }
226
227 // create the ESD output file
228 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
229 if (!file->IsOpen()) {
230 Error("Run", "opening AliESDs.root failed");
231 if (fStopOnError) return kFALSE;
232 }
233
234 // loop over events
235 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
236 Info("Run", "processing event %d", iEvent);
237 AliESD* esd = new AliESD;
238 fRunLoader->GetEvent(iEvent);
239 esd->SetRunNumber(gAlice->GetRunNumber());
240 esd->SetEventNumber(gAlice->GetEvNumber());
241
242 // barrel tracking
243 if (fRunTracking) {
244 if (!RunTracking(esd)) {
245 if (fStopOnError) return kFALSE;
246 }
247 }
248
249 // fill ESD
250 if (!fFillESD.IsNull()) {
251 if (!FillESD(esd, fFillESD)) {
252 if (fStopOnError) return kFALSE;
253 }
254 }
255
256 // combined PID
257 AliESDpid::MakePID(esd);
258
259 // write ESD
260 char name[100];
261 sprintf(name, "ESD%d", iEvent);
262 file->cd();
263 if (!esd->Write(name)) {
264 Error("Run", "writing ESD failed");
265 if (fStopOnError) return kFALSE;
266 }
267 }
268
269 file->Close();
270
271 return kTRUE;
272}
273
274
275//_____________________________________________________________________________
276Bool_t AliReconstruction::RunReconstruction(const TString& detectors)
277{
278// run the reconstruction
279
030b532d 280 TStopwatch stopwatch;
281 stopwatch.Start();
282
596a855f 283 TString detStr = detectors;
284 TObjArray* detArray = gAlice->Detectors();
285 for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
286 AliModule* det = (AliModule*) detArray->At(iDet);
287 if (!det || !det->IsActive()) continue;
288 if (IsSelected(det->GetName(), detStr)) {
289 Info("RunReconstruction", "running reconstruction for %s",
290 det->GetName());
030b532d 291 TStopwatch stopwatchDet;
292 stopwatchDet.Start();
596a855f 293 det->Reconstruct();
030b532d 294 Info("RunReconstruction", "execution time for %s:", det->GetName());
295 stopwatchDet.Print();
596a855f 296 }
297 }
298
299 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
300 Error("RunReconstruction", "the following detectors were not found: %s",
301 detStr.Data());
302 if (fStopOnError) return kFALSE;
303 }
304
030b532d 305 Info("RunReconstruction", "execution time:");
306 stopwatch.Print();
307
596a855f 308 return kTRUE;
309}
310
311//_____________________________________________________________________________
312Bool_t AliReconstruction::RunTracking(AliESD* esd)
313{
314// run the barrel tracking
315
030b532d 316 TStopwatch stopwatch;
317 stopwatch.Start();
318
596a855f 319 // get the primary vertex (from MC for the moment)
320 TArrayF vertex(3);
321 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
322 Double_t vtxPos[3] = {vertex[0], vertex[1], vertex[2]};
323 Double_t vtxCov[6] = {
324 0.005,
325 0.000, 0.005,
326 0.000, 0.000, 0.010
327 };
328 Double_t vtxErr[3] = {vtxCov[0], vtxCov[2], vtxCov[5]}; // diag. elements
329 esd->SetVertex(vtxPos, vtxCov);
330 fITSTracker->SetVertex(vtxPos, vtxErr);
331 fTPCTracker->SetVertex(vtxPos, vtxErr);
332 fTRDTracker->SetVertex(vtxPos, vtxErr);
333
334 // TPC inward fit
335 Info("RunTracking", "TPC inward fit");
336 fTPCLoader->LoadRecPoints("read");
337 TTree* tpcTree = fTPCLoader->TreeR();
338 if (!tpcTree) {
339 Error("RunTracking", "Can't get the TPC cluster tree");
340 return kFALSE;
341 }
342 fTPCTracker->LoadClusters(tpcTree);
343 if (fTPCTracker->Clusters2Tracks(esd) != 0) {
344 Error("RunTracking", "TPC inward track fit failed");
345 return kFALSE;
346 }
347
348 // ITS inward fit
349 Info("RunTracking", "ITS inward fit");
350 fITSLoader->LoadRecPoints("read");
351 TTree* itsTree = fITSLoader->TreeR();
352 if (!itsTree) {
353 Error("RunTracking", "Can't get the ITS cluster tree");
354 return kFALSE;
355 }
356 fITSTracker->LoadClusters(itsTree);
357 if (fITSTracker->Clusters2Tracks(esd) != 0) {
358 Error("RunTracking", "ITS inward track fit failed");
359 return kFALSE;
360 }
361
362 // ITS back propagation
363 Info("RunTracking", "ITS back propagation");
364 if (fITSTracker->PropagateBack(esd) != 0) {
365 Error("RunTracking", "ITS backward propagation failed");
366 return kFALSE;
367 }
368
369 // TPC back propagation
370 Info("RunTracking", "TPC back propagation");
371 if (fTPCTracker->PropagateBack(esd) != 0) {
372 Error("RunTracking", "TPC backward propagation failed");
373 return kFALSE;
374 }
375
376 // TRD back propagation
377 Info("RunTracking", "TRD back propagation");
378 fTRDLoader->LoadRecPoints("read");
379 TTree* trdTree = fTRDLoader->TreeR();
380 if (!trdTree) {
381 Error("RunTracking", "Can't get the TRD cluster tree");
382 return kFALSE;
383 }
384 fTRDTracker->LoadClusters(trdTree);
385 if (fTRDTracker->PropagateBack(esd) != 0) {
386 Error("RunTracking", "TRD backward propagation failed");
387 return kFALSE;
388 }
389
390 // TOF back propagation
391 Info("RunTracking", "TOF back propagation");
392 fTOFLoader->LoadDigits("read");
393 TTree* tofTree = fTOFLoader->TreeD();
394 if (!tofTree) {
395 Error("RunTracking", "Can't get the TOF digits tree");
396 return kFALSE;
397 }
398 fTOFTracker->LoadClusters(tofTree);
399 if (fTOFTracker->PropagateBack(esd) != 0) {
400 Error("RunTracking", "TOF backward propagation failed");
401 return kFALSE;
402 }
403 fTOFTracker->UnloadClusters();
404 fTOFLoader->UnloadDigits();
405
406 // TRD inward refit
407 Info("RunTracking", "TRD inward refit");
408 if (fTRDTracker->RefitInward(esd) != 0) {
409 Error("RunTracking", "TRD inward refit failed");
410 return kFALSE;
411 }
412 fTRDTracker->UnloadClusters();
413 fTRDLoader->UnloadRecPoints();
414
415 // TPC inward refit
416 Info("RunTracking", "TPC inward refit");
417 if (fTPCTracker->RefitInward(esd) != 0) {
418 Error("RunTracking", "TPC inward refit failed");
419 return kFALSE;
420 }
421 fTPCTracker->UnloadClusters();
422 fTPCLoader->UnloadRecPoints();
423
424 // ITS inward refit
425 Info("RunTracking", "ITS inward refit");
426 if (fITSTracker->RefitInward(esd) != 0) {
427 Error("RunTracking", "ITS inward refit failed");
428 return kFALSE;
429 }
430 fITSTracker->UnloadClusters();
431 fITSLoader->UnloadRecPoints();
432
030b532d 433 Info("RunTracking", "execution time:");
434 stopwatch.Print();
435
596a855f 436 return kTRUE;
437}
438
439//_____________________________________________________________________________
440Bool_t AliReconstruction::FillESD(AliESD* esd, const TString& detectors)
441{
442// fill the event summary data
443
030b532d 444 TStopwatch stopwatch;
445 stopwatch.Start();
446
596a855f 447 TString detStr = detectors;
448 TObjArray* detArray = gAlice->Detectors();
449 for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
450 AliModule* det = (AliModule*) detArray->At(iDet);
451 if (!det || !det->IsActive()) continue;
452 if (IsSelected(det->GetName(), detStr)) {
453 Info("FillESD", "filling ESD for %s",
454 det->GetName());
455 det->FillESD(esd);
456 }
457 }
458
459 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
460 Error("FillESD", "the following detectors were not found: %s",
461 detStr.Data());
462 if (fStopOnError) return kFALSE;
463 }
464
030b532d 465 Info("FillESD", "execution time:");
466 stopwatch.Print();
467
596a855f 468 return kTRUE;
469}
470
471
472//_____________________________________________________________________________
473Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
474{
475// check whether detName is contained in detectors
476// if yes, it is removed from detectors
477
478 // check if all detectors are selected
479 if ((detectors.CompareTo("ALL") == 0) ||
480 detectors.BeginsWith("ALL ") ||
481 detectors.EndsWith(" ALL") ||
482 detectors.Contains(" ALL ")) {
483 detectors = "ALL";
484 return kTRUE;
485 }
486
487 // search for the given detector
488 Bool_t result = kFALSE;
489 if ((detectors.CompareTo(detName) == 0) ||
490 detectors.BeginsWith(detName+" ") ||
491 detectors.EndsWith(" "+detName) ||
492 detectors.Contains(" "+detName+" ")) {
493 detectors.ReplaceAll(detName, "");
494 result = kTRUE;
495 }
496
497 // clean up the detectors string
498 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
499 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
500 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
501
502 return result;
503}