9df426ef4dee6914d464886503e062b5e7804032
[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   if (fRunTracking) {
186     // get loaders and trackers
187     fITSLoader = fRunLoader->GetLoader("ITSLoader");
188     if (!fITSLoader) {
189       Error("Run", "no ITS loader found");
190       if (fStopOnError) {CleanUp(); return kFALSE;}
191     }
192     fITSTracker = NULL;
193     if (aliRun->GetDetector("ITS")) {
194       fITSTracker = aliRun->GetDetector("ITS")->CreateTracker();
195     }
196     if (!fITSTracker) {
197       Error("Run", "couldn't create a tracker for ITS");
198       if (fStopOnError) {CleanUp(); return kFALSE;}
199     }
200     
201     fTPCLoader = fRunLoader->GetLoader("TPCLoader");
202     if (!fTPCLoader) {
203       Error("Run", "no TPC loader found");
204       if (fStopOnError) {CleanUp(); return kFALSE;}
205     }
206     fTPCTracker = NULL;
207     if (aliRun->GetDetector("TPC")) {
208       fTPCTracker = aliRun->GetDetector("TPC")->CreateTracker();
209     }
210     if (!fTPCTracker) {
211       Error("Run", "couldn't create a tracker for TPC");
212       if (fStopOnError) {CleanUp(); return kFALSE;}
213     }
214     
215     fTRDLoader = fRunLoader->GetLoader("TRDLoader");
216     if (!fTRDLoader) {
217       Error("Run", "no TRD loader found");
218       if (fStopOnError) {CleanUp(); return kFALSE;}
219     }
220     fTRDTracker = NULL;
221     if (aliRun->GetDetector("TRD")) {
222       fTRDTracker = aliRun->GetDetector("TRD")->CreateTracker();
223     }
224     if (!fTRDTracker) {
225       Error("Run", "couldn't create a tracker for TRD");
226       if (fStopOnError) {CleanUp(); return kFALSE;}
227     }
228     
229     fTOFLoader = fRunLoader->GetLoader("TOFLoader");
230     if (!fTOFLoader) {
231       Error("Run", "no TOF loader found");
232       if (fStopOnError) {CleanUp(); return kFALSE;}
233     }
234     fTOFTracker = NULL;
235     if (aliRun->GetDetector("TOF")) {
236       fTOFTracker = aliRun->GetDetector("TOF")->CreateTracker();
237     }
238     if (!fTOFTracker) {
239       Error("Run", "couldn't create a tracker for TOF");
240       if (fStopOnError) {CleanUp(); return kFALSE;}
241     }
242   }
243   // create the ESD output file
244   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
245   if (!file->IsOpen()) {
246     Error("Run", "opening AliESDs.root failed");
247     if (fStopOnError) {CleanUp(file); return kFALSE;}    
248   }
249
250   // loop over events
251   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
252     Info("Run", "processing event %d", iEvent);
253     AliESD* esd = new AliESD;
254     fRunLoader->GetEvent(iEvent);
255     esd->SetRunNumber(aliRun->GetRunNumber());
256     esd->SetEventNumber(aliRun->GetEvNumber());
257     esd->SetMagneticField(aliRun->Field()->SolenoidField());
258
259     // barrel tracking
260     if (fRunTracking) {
261       if (!RunTracking(esd)) {
262         if (fStopOnError) {CleanUp(file); return kFALSE;}
263       }
264     }
265
266     // fill ESD
267     if (!fFillESD.IsNull()) {
268       if (!FillESD(esd, fFillESD)) {
269         if (fStopOnError) {CleanUp(file); return kFALSE;}
270       }
271     }
272
273     // combined PID
274     AliESDpid::MakePID(esd);
275
276     // write ESD
277     char name[100]; 
278     sprintf(name, "ESD%d", iEvent);
279     file->cd();
280     if (!esd->Write(name)) {
281       Error("Run", "writing ESD failed");
282       if (fStopOnError) {CleanUp(file); return kFALSE;}
283     }
284   }
285
286   CleanUp(file);
287
288   return kTRUE;
289 }
290
291
292 //_____________________________________________________________________________
293 Bool_t AliReconstruction::RunReconstruction(const TString& detectors)
294 {
295 // run the reconstruction
296
297   TStopwatch stopwatch;
298   stopwatch.Start();
299
300   TString detStr = detectors;
301   TObjArray* detArray = fRunLoader->GetAliRun()->Detectors();
302   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
303     AliModule* det = (AliModule*) detArray->At(iDet);
304     if (!det || !det->IsActive()) continue;
305     if (IsSelected(det->GetName(), detStr)) {
306       Info("RunReconstruction", "running reconstruction for %s", 
307            det->GetName());
308       TStopwatch stopwatchDet;
309       stopwatchDet.Start();
310       det->Reconstruct();
311       Info("RunReconstruction", "execution time for %s:", det->GetName());
312       stopwatchDet.Print();
313     }
314   }
315
316   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
317     Error("RunReconstruction", "the following detectors were not found: %s", 
318           detStr.Data());
319     if (fStopOnError) return kFALSE;
320   }
321
322   Info("RunReconstruction", "execution time:");
323   stopwatch.Print();
324
325   return kTRUE;
326 }
327
328 //_____________________________________________________________________________
329 Bool_t AliReconstruction::RunTracking(AliESD* esd)
330 {
331 // run the barrel tracking
332
333   TStopwatch stopwatch;
334   stopwatch.Start();
335
336   // get the primary vertex (from MC for the moment)
337   TArrayF vertex(3);     
338   fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
339   Double_t vtxPos[3] = {vertex[0], vertex[1], vertex[2]};
340   Double_t vtxCov[6] = {
341     0.005,
342     0.000, 0.005,
343     0.000, 0.000, 0.010
344   };
345   Double_t vtxErr[3] = {vtxCov[0], vtxCov[2], vtxCov[5]}; // diag. elements
346   esd->SetVertex(vtxPos, vtxCov);
347   fITSTracker->SetVertex(vtxPos, vtxErr);
348   fTPCTracker->SetVertex(vtxPos, vtxErr);
349   fTRDTracker->SetVertex(vtxPos, vtxErr);
350
351   // TPC tracking
352   Info("RunTracking", "TPC tracking");
353   fTPCLoader->LoadRecPoints("read");
354   TTree* tpcTree = fTPCLoader->TreeR();
355   if (!tpcTree) {
356     Error("RunTracking", "Can't get the TPC cluster tree");
357     return kFALSE;
358   }     
359   fTPCTracker->LoadClusters(tpcTree);
360   if (fTPCTracker->Clusters2Tracks(esd) != 0) {
361     Error("RunTracking", "TPC Clusters2Tracks failed");
362     return kFALSE;
363   }
364
365   fRunLoader->GetAliRun()->GetDetector("TPC")->FillESD(esd); // preliminary PID
366   AliESDpid::MakePID(esd);                  // for the ITS tracker
367
368   // ITS tracking
369   Info("RunTracking", "ITS tracking");
370   fITSLoader->LoadRecPoints("read");
371   TTree* itsTree = fITSLoader->TreeR();
372   if (!itsTree) {
373     Error("RunTracking", "Can't get the ITS cluster tree");
374     return kFALSE;
375   }     
376   fITSTracker->LoadClusters(itsTree);
377   if (fITSTracker->Clusters2Tracks(esd) != 0) {
378     Error("RunTracking", "ITS Clusters2Tracks failed");
379     return kFALSE;
380   }
381
382   // ITS back propagation
383   Info("RunTracking", "ITS back propagation");
384   if (fITSTracker->PropagateBack(esd) != 0) {
385     Error("RunTracking", "ITS backward propagation failed");
386     return kFALSE;
387   }
388
389   // TPC back propagation
390   Info("RunTracking", "TPC back propagation");
391   if (fTPCTracker->PropagateBack(esd) != 0) {
392     Error("RunTracking", "TPC backward propagation failed");
393     return kFALSE;
394   }
395
396   // TRD back propagation
397   Info("RunTracking", "TRD back propagation");
398   fTRDLoader->LoadRecPoints("read");
399   TTree* trdTree = fTRDLoader->TreeR();
400   if (!trdTree) {
401     Error("RunTracking", "Can't get the TRD cluster tree");
402     return kFALSE;
403   }     
404   fTRDTracker->LoadClusters(trdTree);
405   if (fTRDTracker->PropagateBack(esd) != 0) {
406     Error("RunTracking", "TRD backward propagation failed");
407     return kFALSE;
408   }
409
410   // TOF back propagation
411   Info("RunTracking", "TOF back propagation");
412   fTOFLoader->LoadDigits("read");
413   TTree* tofTree = fTOFLoader->TreeD();
414   if (!tofTree) {
415     Error("RunTracking", "Can't get the TOF digits tree");
416     return kFALSE;
417   }     
418   fTOFTracker->LoadClusters(tofTree);
419   if (fTOFTracker->PropagateBack(esd) != 0) {
420     Error("RunTracking", "TOF backward propagation failed");
421     return kFALSE;
422   }
423   fTOFTracker->UnloadClusters();
424   fTOFLoader->UnloadDigits();
425
426   // TRD inward refit
427   Info("RunTracking", "TRD inward refit");
428   if (fTRDTracker->RefitInward(esd) != 0) {
429     Error("RunTracking", "TRD inward refit failed");
430     return kFALSE;
431   }
432   fTRDTracker->UnloadClusters();
433   fTRDLoader->UnloadRecPoints();
434     
435   // TPC inward refit
436   Info("RunTracking", "TPC inward refit");
437   if (fTPCTracker->RefitInward(esd) != 0) {
438     Error("RunTracking", "TPC inward refit failed");
439     return kFALSE;
440   }
441   fTPCTracker->UnloadClusters();
442   fTPCLoader->UnloadRecPoints();
443     
444   // ITS inward refit
445   Info("RunTracking", "ITS inward refit");
446   if (fITSTracker->RefitInward(esd) != 0) {
447     Error("RunTracking", "ITS inward refit failed");
448     return kFALSE;
449   }
450   fITSTracker->UnloadClusters();
451   fITSLoader->UnloadRecPoints();
452
453   Info("RunTracking", "execution time:");
454   stopwatch.Print();
455
456   return kTRUE;
457 }
458
459 //_____________________________________________________________________________
460 Bool_t AliReconstruction::FillESD(AliESD* esd, const TString& detectors)
461 {
462 // fill the event summary data
463
464   TStopwatch stopwatch;
465   stopwatch.Start();
466
467   TString detStr = detectors;
468   TObjArray* detArray = fRunLoader->GetAliRun()->Detectors();
469   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
470     AliModule* det = (AliModule*) detArray->At(iDet);
471     if (!det || !det->IsActive()) continue;
472     if (IsSelected(det->GetName(), detStr)) {
473       Info("FillESD", "filling ESD for %s", 
474            det->GetName());
475       det->FillESD(esd);
476     }
477   }
478
479   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
480     Error("FillESD", "the following detectors were not found: %s", 
481           detStr.Data());
482     if (fStopOnError) return kFALSE;
483   }
484
485   Info("FillESD", "execution time:");
486   stopwatch.Print();
487
488   return kTRUE;
489 }
490
491
492 //_____________________________________________________________________________
493 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
494 {
495 // check whether detName is contained in detectors
496 // if yes, it is removed from detectors
497
498   // check if all detectors are selected
499   if ((detectors.CompareTo("ALL") == 0) ||
500       detectors.BeginsWith("ALL ") ||
501       detectors.EndsWith(" ALL") ||
502       detectors.Contains(" ALL ")) {
503     detectors = "ALL";
504     return kTRUE;
505   }
506
507   // search for the given detector
508   Bool_t result = kFALSE;
509   if ((detectors.CompareTo(detName) == 0) ||
510       detectors.BeginsWith(detName+" ") ||
511       detectors.EndsWith(" "+detName) ||
512       detectors.Contains(" "+detName+" ")) {
513     detectors.ReplaceAll(detName, "");
514     result = kTRUE;
515   }
516
517   // clean up the detectors string
518   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
519   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
520   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
521
522   return result;
523 }
524
525 //_____________________________________________________________________________
526 void AliReconstruction::CleanUp(TFile* file)
527 {
528 // delete trackers and the run loader and close and delete the file
529
530   delete fITSTracker;
531   fITSTracker = NULL;
532   delete fTPCTracker;
533   fTPCTracker = NULL;
534   delete fTRDTracker;
535   fTRDTracker = NULL;
536   delete fTOFTracker;
537   fTOFTracker = NULL;
538
539   delete fRunLoader;
540   fRunLoader = NULL;
541
542   if (file) {
543     file->Close();
544     delete file;
545   }
546 }