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