d7c1c878a3655fddbd1c815dbb959818b9c429ad
[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 // For debug purposes the method SetCheckPointLevel can be used. If the      //
58 // argument is greater than 0, files with ESD events will be written after   //
59 // selected steps of the reconstruction for each event:                      //
60 //   level 1: after tracking and after filling of ESD (final)                //
61 //   level 2: in addition after each tracking step                           //
62 //   level 3: in addition after the filling of ESD for each detector         //
63 // If a final check point file exists for an event, this event will be       //
64 // skipped in the reconstruction. The tracking and the filling of ESD for    //
65 // a detector will be skipped as well, if the corresponding check point      //
66 // file exists. The ESD event will then be loaded from the file instead.     //
67 //                                                                           //
68 ///////////////////////////////////////////////////////////////////////////////
69
70
71 #include "AliReconstruction.h"
72 #include "AliRunLoader.h"
73 #include "AliRun.h"
74 #include "AliModule.h"
75 #include "AliDetector.h"
76 #include "AliTracker.h"
77 #include "AliESD.h"
78 #include "AliHeader.h"
79 #include "AliGenEventHeader.h"
80 #include "AliESDpid.h"
81 #include "AliMagF.h"
82 #include <TArrayF.h>
83 #include <TSystem.h>
84 #include <TROOT.h>
85
86
87 ClassImp(AliReconstruction)
88
89
90 //_____________________________________________________________________________
91 AliReconstruction::AliReconstruction(const char* gAliceFilename,
92                                      const char* name, const char* title) :
93   TNamed(name, title),
94
95   fRunReconstruction("ALL"),
96   fRunTracking(kTRUE),
97   fFillESD("ALL"),
98   fGAliceFileName(gAliceFilename),
99   fStopOnError(kFALSE),
100   fCheckPointLevel(0),
101
102   fRunLoader(NULL),
103   fITSLoader(NULL),
104   fITSTracker(NULL),
105   fTPCLoader(NULL),
106   fTPCTracker(NULL),
107   fTRDLoader(NULL),
108   fTRDTracker(NULL),
109   fTOFLoader(NULL),
110   fTOFTracker(NULL)
111 {
112 // create reconstruction object with default parameters
113
114 }
115
116 //_____________________________________________________________________________
117 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
118   TNamed(rec),
119
120   fRunReconstruction(rec.fRunReconstruction),
121   fRunTracking(rec.fRunTracking),
122   fFillESD(rec.fFillESD),
123   fGAliceFileName(rec.fGAliceFileName),
124   fStopOnError(rec.fStopOnError),
125   fCheckPointLevel(0),
126
127   fRunLoader(NULL),
128   fITSLoader(NULL),
129   fITSTracker(NULL),
130   fTPCLoader(NULL),
131   fTPCTracker(NULL),
132   fTRDLoader(NULL),
133   fTRDTracker(NULL),
134   fTOFLoader(NULL),
135   fTOFTracker(NULL)
136 {
137 // copy constructor
138
139 }
140
141 //_____________________________________________________________________________
142 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
143 {
144 // assignment operator
145
146   this->~AliReconstruction();
147   new(this) AliReconstruction(rec);
148   return *this;
149 }
150
151 //_____________________________________________________________________________
152 AliReconstruction::~AliReconstruction()
153 {
154 // clean up
155
156   CleanUp();
157 }
158
159
160 //_____________________________________________________________________________
161 void AliReconstruction::SetGAliceFile(const char* fileName)
162 {
163 // set the name of the galice file
164
165   fGAliceFileName = fileName;
166 }
167
168
169 //_____________________________________________________________________________
170 Bool_t AliReconstruction::Run()
171 {
172 // run the reconstruction
173
174   // open the run loader
175   fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
176   if (!fRunLoader) {
177     Error("Run", "no run loader found in file %s", 
178           fGAliceFileName.Data());
179     CleanUp();
180     return kFALSE;
181   }
182   fRunLoader->LoadgAlice();
183   AliRun* aliRun = fRunLoader->GetAliRun();
184   if (!aliRun) {
185     Error("Run", "no gAlice object found in file %s", 
186           fGAliceFileName.Data());
187     CleanUp();
188     return kFALSE;
189   }
190   gAlice = aliRun;
191
192   // local reconstruction
193   if (!fRunReconstruction.IsNull()) {
194     if (!RunReconstruction(fRunReconstruction)) {
195       if (fStopOnError) {CleanUp(); return kFALSE;}
196     }
197   }
198   if (!fRunTracking && fFillESD.IsNull()) return kTRUE;
199
200   // get loaders and trackers
201   if (fRunTracking && !CreateTrackers()) {
202     if (fStopOnError) {
203       CleanUp(); 
204       return kFALSE;
205     }      
206   }
207
208   // create the ESD output file
209   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
210   if (!file->IsOpen()) {
211     Error("Run", "opening AliESDs.root failed");
212     if (fStopOnError) {CleanUp(file); return kFALSE;}    
213   }
214
215   // loop over events
216   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
217     Info("Run", "processing event %d", iEvent);
218     fRunLoader->GetEvent(iEvent);
219
220     char fileName[256];
221     sprintf(fileName, "ESD_%d.%d_final.root", 
222             aliRun->GetRunNumber(), aliRun->GetEvNumber());
223     if (!gSystem->AccessPathName(fileName)) continue;
224
225     AliESD* esd = new AliESD;
226     esd->SetRunNumber(aliRun->GetRunNumber());
227     esd->SetEventNumber(aliRun->GetEvNumber());
228     esd->SetMagneticField(aliRun->Field()->SolenoidField());
229
230     // barrel tracking
231     if (fRunTracking) {
232       if (!ReadESD(esd, "tracking")) {
233         if (!RunTracking(esd)) {
234           if (fStopOnError) {CleanUp(file); return kFALSE;}
235         }
236         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
237       }
238     }
239
240     // fill ESD
241     if (!fFillESD.IsNull()) {
242       if (!FillESD(esd, fFillESD)) {
243         if (fStopOnError) {CleanUp(file); return kFALSE;}
244       }
245     }
246
247     // combined PID
248     AliESDpid::MakePID(esd);
249     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
250
251     // write ESD
252     char name[100]; 
253     sprintf(name, "ESD%d", iEvent);
254     file->cd();
255     if (!esd->Write(name)) {
256       Error("Run", "writing ESD failed");
257       if (fStopOnError) {CleanUp(file); return kFALSE;}
258     }
259     file->Flush();
260
261     if (fCheckPointLevel > 0) WriteESD(esd, "final");
262     delete esd;
263   }
264
265   CleanUp(file);
266
267   return kTRUE;
268 }
269
270
271 //_____________________________________________________________________________
272 Bool_t AliReconstruction::RunReconstruction(const TString& detectors)
273 {
274 // run the reconstruction
275
276   TStopwatch stopwatch;
277   stopwatch.Start();
278
279   TString detStr = detectors;
280   TObjArray* detArray = fRunLoader->GetAliRun()->Detectors();
281   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
282     AliModule* det = (AliModule*) detArray->At(iDet);
283     if (!det || !det->IsActive()) continue;
284     if (IsSelected(det->GetName(), detStr)) {
285       Info("RunReconstruction", "running reconstruction for %s", 
286            det->GetName());
287       TStopwatch stopwatchDet;
288       stopwatchDet.Start();
289       det->Reconstruct();
290       Info("RunReconstruction", "execution time for %s:", det->GetName());
291       stopwatchDet.Print();
292     }
293   }
294
295   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
296     Error("RunReconstruction", "the following detectors were not found: %s", 
297           detStr.Data());
298     if (fStopOnError) return kFALSE;
299   }
300
301   Info("RunReconstruction", "execution time:");
302   stopwatch.Print();
303
304   return kTRUE;
305 }
306
307 //_____________________________________________________________________________
308 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
309 {
310 // run the barrel tracking
311
312   TStopwatch stopwatch;
313   stopwatch.Start();
314
315   // get the primary vertex (from MC for the moment)
316   TArrayF vertex(3);     
317   fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
318   Double_t vtxPos[3] = {vertex[0], vertex[1], vertex[2]};
319   Double_t vtxCov[6] = {
320     0.005,
321     0.000, 0.005,
322     0.000, 0.000, 0.010
323   };
324   Double_t vtxErr[3] = {vtxCov[0], vtxCov[2], vtxCov[5]}; // diag. elements
325   esd->SetVertex(vtxPos, vtxCov);
326   if (fITSTracker) fITSTracker->SetVertex(vtxPos, vtxErr);
327   if (fTPCTracker) fTPCTracker->SetVertex(vtxPos, vtxErr);
328   if (fTRDTracker) fTRDTracker->SetVertex(vtxPos, vtxErr);
329   if (fCheckPointLevel > 1) WriteESD(esd, "vertex");
330
331   if (!fTPCTracker) {
332     Error("RunTracking", "no TPC tracker");
333     return kFALSE;
334   }
335
336   // TPC tracking
337   Info("RunTracking", "TPC tracking");
338   fTPCLoader->LoadRecPoints("read");
339   TTree* tpcTree = fTPCLoader->TreeR();
340   if (!tpcTree) {
341     Error("RunTracking", "Can't get the TPC cluster tree");
342     return kFALSE;
343   }
344   fTPCTracker->LoadClusters(tpcTree);
345   if (fTPCTracker->Clusters2Tracks(esd) != 0) {
346     Error("RunTracking", "TPC Clusters2Tracks failed");
347     return kFALSE;
348   }
349   if (fCheckPointLevel > 1) WriteESD(esd, "TPC.tracking");
350
351   if (!fITSTracker) {
352     Warning("RunTracking", "no ITS tracker");
353   } else {
354
355     fRunLoader->GetAliRun()->GetDetector("TPC")->FillESD(esd); // preliminary
356     AliESDpid::MakePID(esd);                  // PID for the ITS tracker
357
358     // ITS tracking
359     Info("RunTracking", "ITS tracking");
360     fITSLoader->LoadRecPoints("read");
361     TTree* itsTree = fITSLoader->TreeR();
362     if (!itsTree) {
363       Error("RunTracking", "Can't get the ITS cluster tree");
364       return kFALSE;
365     }
366     fITSTracker->LoadClusters(itsTree);
367     if (fITSTracker->Clusters2Tracks(esd) != 0) {
368       Error("RunTracking", "ITS Clusters2Tracks failed");
369       return kFALSE;
370     }
371     if (fCheckPointLevel > 1) WriteESD(esd, "ITS.tracking");
372
373     if (!fTRDTracker) {
374       Warning("RunTracking", "no TRD tracker");
375     } else {
376       // ITS back propagation
377       Info("RunTracking", "ITS back propagation");
378       if (fITSTracker->PropagateBack(esd) != 0) {
379         Error("RunTracking", "ITS backward propagation failed");
380         return kFALSE;
381       }
382       if (fCheckPointLevel > 1) WriteESD(esd, "ITS.back");
383
384       // TPC back propagation
385       Info("RunTracking", "TPC back propagation");
386       if (fTPCTracker->PropagateBack(esd) != 0) {
387         Error("RunTracking", "TPC backward propagation failed");
388         return kFALSE;
389       }
390       if (fCheckPointLevel > 1) WriteESD(esd, "TPC.back");
391
392       // TRD back propagation
393       Info("RunTracking", "TRD back propagation");
394       fTRDLoader->LoadRecPoints("read");
395       TTree* trdTree = fTRDLoader->TreeR();
396       if (!trdTree) {
397         Error("RunTracking", "Can't get the TRD cluster tree");
398         return kFALSE;
399       }
400       fTRDTracker->LoadClusters(trdTree);
401       if (fTRDTracker->PropagateBack(esd) != 0) {
402         Error("RunTracking", "TRD backward propagation failed");
403         return kFALSE;
404       }
405       if (fCheckPointLevel > 1) WriteESD(esd, "TRD.back");
406
407       if (!fTOFTracker) {
408         Warning("RunTracking", "no TOF tracker");
409       } else {
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         if (fCheckPointLevel > 1) WriteESD(esd, "TOF.back");
424         fTOFTracker->UnloadClusters();
425         fTOFLoader->UnloadDigits();
426       }
427
428       // TRD inward refit
429       Info("RunTracking", "TRD inward refit");
430       if (fTRDTracker->RefitInward(esd) != 0) {
431         Error("RunTracking", "TRD inward refit failed");
432         return kFALSE;
433       }
434       if (fCheckPointLevel > 1) WriteESD(esd, "TRD.refit");
435       fTRDTracker->UnloadClusters();
436       fTRDLoader->UnloadRecPoints();
437     
438       // TPC inward refit
439       Info("RunTracking", "TPC inward refit");
440       if (fTPCTracker->RefitInward(esd) != 0) {
441         Error("RunTracking", "TPC inward refit failed");
442         return kFALSE;
443       }
444       if (fCheckPointLevel > 1) WriteESD(esd, "TPC.refit");
445     
446       // ITS inward refit
447       Info("RunTracking", "ITS inward refit");
448       if (fITSTracker->RefitInward(esd) != 0) {
449         Error("RunTracking", "ITS inward refit failed");
450         return kFALSE;
451       }
452       if (fCheckPointLevel > 1) WriteESD(esd, "ITS.refit");
453
454     }  // if TRD tracker
455     fITSTracker->UnloadClusters();
456     fITSLoader->UnloadRecPoints();
457
458   }  // if ITS tracker
459   fTPCTracker->UnloadClusters();
460   fTPCLoader->UnloadRecPoints();
461
462   Info("RunTracking", "execution time:");
463   stopwatch.Print();
464
465   return kTRUE;
466 }
467
468 //_____________________________________________________________________________
469 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
470 {
471 // fill the event summary data
472
473   TStopwatch stopwatch;
474   stopwatch.Start();
475
476   TString detStr = detectors;
477   TObjArray* detArray = fRunLoader->GetAliRun()->Detectors();
478   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
479     AliModule* det = (AliModule*) detArray->At(iDet);
480     if (!det || !det->IsActive()) continue;
481     if (IsSelected(det->GetName(), detStr)) {
482       if (!ReadESD(esd, det->GetName())) {
483         Info("FillESD", "filling ESD for %s", 
484              det->GetName());
485         det->FillESD(esd);
486         if (fCheckPointLevel > 2) WriteESD(esd, det->GetName());
487       }
488     }
489   }
490
491   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
492     Error("FillESD", "the following detectors were not found: %s", 
493           detStr.Data());
494     if (fStopOnError) return kFALSE;
495   }
496
497   Info("FillESD", "execution time:");
498   stopwatch.Print();
499
500   return kTRUE;
501 }
502
503
504 //_____________________________________________________________________________
505 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
506 {
507 // check whether detName is contained in detectors
508 // if yes, it is removed from detectors
509
510   // check if all detectors are selected
511   if ((detectors.CompareTo("ALL") == 0) ||
512       detectors.BeginsWith("ALL ") ||
513       detectors.EndsWith(" ALL") ||
514       detectors.Contains(" ALL ")) {
515     detectors = "ALL";
516     return kTRUE;
517   }
518
519   // search for the given detector
520   Bool_t result = kFALSE;
521   if ((detectors.CompareTo(detName) == 0) ||
522       detectors.BeginsWith(detName+" ") ||
523       detectors.EndsWith(" "+detName) ||
524       detectors.Contains(" "+detName+" ")) {
525     detectors.ReplaceAll(detName, "");
526     result = kTRUE;
527   }
528
529   // clean up the detectors string
530   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
531   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
532   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
533
534   return result;
535 }
536
537 //_____________________________________________________________________________
538 Bool_t AliReconstruction::CreateTrackers()
539 {
540 // get the loaders and create the trackers
541
542   AliRun* aliRun = fRunLoader->GetAliRun();
543
544   fITSTracker = NULL;
545   fITSLoader = fRunLoader->GetLoader("ITSLoader");
546   if (!fITSLoader) {
547     Warning("CreateTrackers", "no ITS loader found");
548     if (fStopOnError) return kFALSE;
549   } else {
550     if (aliRun->GetDetector("ITS")) {
551       fITSTracker = aliRun->GetDetector("ITS")->CreateTracker();
552     }
553     if (!fITSTracker) {
554       Warning("CreateTrackers", "couldn't create a tracker for ITS");
555       if (fStopOnError) return kFALSE;
556     }
557   }
558     
559   fTPCTracker = NULL;
560   fTPCLoader = fRunLoader->GetLoader("TPCLoader");
561   if (!fTPCLoader) {
562     Error("CreateTrackers", "no TPC loader found");
563     if (fStopOnError) return kFALSE;
564   } else {
565     if (aliRun->GetDetector("TPC")) {
566       fTPCTracker = aliRun->GetDetector("TPC")->CreateTracker();
567     }
568     if (!fTPCTracker) {
569       Error("CreateTrackers", "couldn't create a tracker for TPC");
570       if (fStopOnError) return kFALSE;
571     }
572   }
573     
574   fTRDTracker = NULL;
575   fTRDLoader = fRunLoader->GetLoader("TRDLoader");
576   if (!fTRDLoader) {
577     Warning("CreateTrackers", "no TRD loader found");
578     if (fStopOnError) return kFALSE;
579   } else {
580     if (aliRun->GetDetector("TRD")) {
581       fTRDTracker = aliRun->GetDetector("TRD")->CreateTracker();
582     }
583     if (!fTRDTracker) {
584       Warning("CreateTrackers", "couldn't create a tracker for TRD");
585       if (fStopOnError) return kFALSE;
586     }
587   }
588     
589   fTOFTracker = NULL;
590   fTOFLoader = fRunLoader->GetLoader("TOFLoader");
591   if (!fTOFLoader) {
592     Warning("CreateTrackers", "no TOF loader found");
593     if (fStopOnError) return kFALSE;
594   } else {
595     if (aliRun->GetDetector("TOF")) {
596       fTOFTracker = aliRun->GetDetector("TOF")->CreateTracker();
597     }
598     if (!fTOFTracker) {
599       Warning("CreateTrackers", "couldn't create a tracker for TOF");
600       if (fStopOnError) return kFALSE;
601     }
602   }
603
604   return kTRUE;
605 }
606
607 //_____________________________________________________________________________
608 void AliReconstruction::CleanUp(TFile* file)
609 {
610 // delete trackers and the run loader and close and delete the file
611
612   delete fITSTracker;
613   fITSTracker = NULL;
614   delete fTPCTracker;
615   fTPCTracker = NULL;
616   delete fTRDTracker;
617   fTRDTracker = NULL;
618   delete fTOFTracker;
619   fTOFTracker = NULL;
620
621   delete fRunLoader;
622   fRunLoader = NULL;
623
624   if (file) {
625     file->Close();
626     delete file;
627   }
628 }
629
630
631 //_____________________________________________________________________________
632 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
633 {
634 // read the ESD event from a file
635
636   if (!esd) return kFALSE;
637   char fileName[256];
638   sprintf(fileName, "ESD_%d.%d_%s.root", 
639           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
640   if (gSystem->AccessPathName(fileName)) return kFALSE;
641
642   Info("ReadESD", "reading ESD from file %s", fileName);
643   TFile* file = TFile::Open(fileName);
644   if (!file || !file->IsOpen()) {
645     Error("ReadESD", "opening %s failed", fileName);
646     delete file;
647     return kFALSE;
648   }
649
650   gROOT->cd();
651   delete esd;
652   esd = (AliESD*) file->Get("ESD");
653   file->Close();
654   delete file;
655   return kTRUE;
656 }
657
658 //_____________________________________________________________________________
659 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
660 {
661 // write the ESD event to a file
662
663   if (!esd) return;
664   char fileName[256];
665   sprintf(fileName, "ESD_%d.%d_%s.root", 
666           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
667
668   Info("WriteESD", "writing ESD to file %s", fileName);
669   TFile* file = TFile::Open(fileName, "recreate");
670   if (!file || !file->IsOpen()) {
671     Error("WriteESD", "opening %s failed", fileName);
672   } else {
673     esd->Write("ESD");
674     file->Close();
675   }
676   delete file;
677 }