Further integration of the HLT reconstruction (C.Cheshkov)
[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 //                                                                           //
30 // If the input to the reconstruction are not simulated digits but raw data, //
31 // this can be specified by an argument of the Run method or by the method   //
32 //                                                                           //
33 //   rec.SetInput("...");                                                    //
34 //                                                                           //
35 // The input formats and the corresponding argument are:                     //
36 // - DDL raw data files: directory name, ends with "/"                       //
37 // - raw data root file: root file name, extension ".root"                   //
38 // - raw data DATE file: DATE file name, any other non-empty string          //
39 // - MC root files     : empty string, default                               //
40 //                                                                           //
41 // By default all events are reconstructed. The reconstruction can be        //
42 // limited to a range of events by giving the index of the first and the     //
43 // last event as an argument to the Run method or by calling                 //
44 //                                                                           //
45 //   rec.SetEventRange(..., ...);                                            //
46 //                                                                           //
47 // The index -1 (default) can be used for the last event to indicate no      //
48 // upper limit of the event range.                                           //
49 //                                                                           //
50 // The name of the galice file can be changed from the default               //
51 // "galice.root" by passing it as argument to the AliReconstruction          //
52 // constructor or by                                                         //
53 //                                                                           //
54 //   rec.SetGAliceFile("...");                                               //
55 //                                                                           //
56 // The local reconstruction can be switched on or off for individual         //
57 // detectors by                                                              //
58 //                                                                           //
59 //   rec.SetRunLocalReconstruction("...");                                   //
60 //                                                                           //
61 // The argument is a (case sensitive) string with the names of the           //
62 // detectors separated by a space. The special string "ALL" selects all      //
63 // available detectors. This is the default.                                 //
64 //                                                                           //
65 // The reconstruction of the primary vertex position can be switched off by  //
66 //                                                                           //
67 //   rec.SetRunVertexFinder(kFALSE);                                         //
68 //                                                                           //
69 // The tracking and the creation of ESD tracks can be switched on for        //
70 // selected detectors by                                                     //
71 //                                                                           //
72 //   rec.SetRunTracking("...");                                              //
73 //                                                                           //
74 // The filling of additional ESD information can be steered by               //
75 //                                                                           //
76 //   rec.SetFillESD("...");                                                  //
77 //                                                                           //
78 // Again, for both methods the string specifies the list of detectors.       //
79 // The default is "ALL".                                                     //
80 //                                                                           //
81 // The call of the shortcut method                                           //
82 //                                                                           //
83 //   rec.SetRunReconstruction("...");                                        //
84 //                                                                           //
85 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
86 // SetFillESD with the same detector selecting string as argument.           //
87 //                                                                           //
88 // The reconstruction requires digits or raw data as input. For the creation //
89 // of digits and raw data have a look at the class AliSimulation.            //
90 //                                                                           //
91 // For debug purposes the method SetCheckPointLevel can be used. If the      //
92 // argument is greater than 0, files with ESD events will be written after   //
93 // selected steps of the reconstruction for each event:                      //
94 //   level 1: after tracking and after filling of ESD (final)                //
95 //   level 2: in addition after each tracking step                           //
96 //   level 3: in addition after the filling of ESD for each detector         //
97 // If a final check point file exists for an event, this event will be       //
98 // skipped in the reconstruction. The tracking and the filling of ESD for    //
99 // a detector will be skipped as well, if the corresponding check point      //
100 // file exists. The ESD event will then be loaded from the file instead.     //
101 //                                                                           //
102 ///////////////////////////////////////////////////////////////////////////////
103
104 #include <TArrayF.h>
105 #include <TFile.h>
106 #include <TSystem.h>
107 #include <TROOT.h>
108 #include <TPluginManager.h>
109 #include <TStopwatch.h>
110
111 #include "AliReconstruction.h"
112 #include "AliReconstructor.h"
113 #include "AliLog.h"
114 #include "AliRunLoader.h"
115 #include "AliRun.h"
116 #include "AliRawReaderFile.h"
117 #include "AliRawReaderDate.h"
118 #include "AliRawReaderRoot.h"
119 #include "AliTracker.h"
120 #include "AliESD.h"
121 #include "AliESDVertex.h"
122 #include "AliVertexer.h"
123 #include "AliHeader.h"
124 #include "AliGenEventHeader.h"
125 #include "AliPID.h"
126 #include "AliESDpid.h"
127 #include "AliMagF.h"
128
129 ClassImp(AliReconstruction)
130
131
132 //_____________________________________________________________________________
133 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
134
135 //_____________________________________________________________________________
136 AliReconstruction::AliReconstruction(const char* gAliceFilename,
137                                      const char* name, const char* title) :
138   TNamed(name, title),
139
140   fRunLocalReconstruction("ALL"),
141   fRunVertexFinder(kTRUE),
142   fRunHLTTracking(kFALSE),
143   fRunTracking("ALL"),
144   fFillESD("ALL"),
145   fGAliceFileName(gAliceFilename),
146   fInput(""),
147   fFirstEvent(0),
148   fLastEvent(-1),
149   fStopOnError(kFALSE),
150   fCheckPointLevel(0),
151   fOptions(),
152
153   fRunLoader(NULL),
154   fRawReader(NULL),
155
156   fVertexer(NULL)
157 {
158 // create reconstruction object with default parameters
159   
160   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
161     fReconstructor[iDet] = NULL;
162     fLoader[iDet] = NULL;
163     fTracker[iDet] = NULL;
164   }
165   //  AliPID::Init();
166 }
167
168 //_____________________________________________________________________________
169 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
170   TNamed(rec),
171
172   fRunLocalReconstruction(rec.fRunLocalReconstruction),
173   fRunVertexFinder(rec.fRunVertexFinder),
174   fRunHLTTracking(rec.fRunHLTTracking),
175   fRunTracking(rec.fRunTracking),
176   fFillESD(rec.fFillESD),
177   fGAliceFileName(rec.fGAliceFileName),
178   fInput(rec.fInput),
179   fFirstEvent(rec.fFirstEvent),
180   fLastEvent(rec.fLastEvent),
181   fStopOnError(rec.fStopOnError),
182   fCheckPointLevel(0),
183   fOptions(),
184
185   fRunLoader(NULL),
186   fRawReader(NULL),
187
188   fVertexer(NULL)
189 {
190 // copy constructor
191
192   for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
193     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
194   }
195   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
196     fReconstructor[iDet] = NULL;
197     fLoader[iDet] = NULL;
198     fTracker[iDet] = NULL;
199   }
200 }
201
202 //_____________________________________________________________________________
203 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
204 {
205 // assignment operator
206
207   this->~AliReconstruction();
208   new(this) AliReconstruction(rec);
209   return *this;
210 }
211
212 //_____________________________________________________________________________
213 AliReconstruction::~AliReconstruction()
214 {
215 // clean up
216
217   CleanUp();
218   fOptions.Delete();
219 }
220
221
222 //_____________________________________________________________________________
223 void AliReconstruction::SetGAliceFile(const char* fileName)
224 {
225 // set the name of the galice file
226
227   fGAliceFileName = fileName;
228 }
229
230 //_____________________________________________________________________________
231 void AliReconstruction::SetOption(const char* detector, const char* option)
232 {
233 // set options for the reconstruction of a detector
234
235   TObject* obj = fOptions.FindObject(detector);
236   if (obj) fOptions.Remove(obj);
237   fOptions.Add(new TNamed(detector, option));
238 }
239
240
241 //_____________________________________________________________________________
242 Bool_t AliReconstruction::Run(const char* input,
243                               Int_t firstEvent, Int_t lastEvent)
244 {
245 // run the reconstruction
246
247   // set the input
248   if (!input) input = fInput.Data();
249   TString fileName(input);
250   if (fileName.EndsWith("/")) {
251     fRawReader = new AliRawReaderFile(fileName);
252   } else if (fileName.EndsWith(".root")) {
253     fRawReader = new AliRawReaderRoot(fileName);
254   } else if (!fileName.IsNull()) {
255     fRawReader = new AliRawReaderDate(fileName);
256     fRawReader->SelectEvents(7);
257   }
258
259   // get the run loader
260   if (!InitRunLoader()) return kFALSE;
261
262   // local reconstruction
263   if (!fRunLocalReconstruction.IsNull()) {
264     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
265       if (fStopOnError) {CleanUp(); return kFALSE;}
266     }
267   }
268 //  if (!fRunVertexFinder && fRunTracking.IsNull() && 
269 //      fFillESD.IsNull()) return kTRUE;
270
271   // get vertexer
272   if (fRunVertexFinder && !CreateVertexer()) {
273     if (fStopOnError) {
274       CleanUp(); 
275       return kFALSE;
276     }
277   }
278
279   // get trackers
280   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
281     if (fStopOnError) {
282       CleanUp(); 
283       return kFALSE;
284     }      
285   }
286
287   // get the possibly already existing ESD file and tree
288   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
289   TFile* fileOld = NULL;
290   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
291   if (!gSystem->AccessPathName("AliESDs.root")){
292     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
293     fileOld = TFile::Open("AliESDs.old.root");
294     if (fileOld && fileOld->IsOpen()) {
295       treeOld = (TTree*) fileOld->Get("esdTree");
296       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
297       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
298       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
299     }
300   }
301
302   // create the ESD output file and tree
303   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
304   if (!file->IsOpen()) {
305     AliError("opening AliESDs.root failed");
306     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}    
307   }
308   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
309   tree->Branch("ESD", "AliESD", &esd);
310   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
311   hlttree->Branch("ESD", "AliESD", &hltesd);
312   delete esd; delete hltesd;
313   esd = NULL; hltesd = NULL;
314   gROOT->cd();
315
316   // loop over events
317   if (fRawReader) fRawReader->RewindEvents();
318   
319   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
320     if (fRawReader) fRawReader->NextEvent();
321     if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
322       // copy old ESD to the new one
323       if (treeOld) {
324         treeOld->SetBranchAddress("ESD", &esd);
325         treeOld->GetEntry(iEvent);
326       }
327       tree->Fill();
328       if (hlttreeOld) {
329         hlttreeOld->SetBranchAddress("ESD", &hltesd);
330         hlttreeOld->GetEntry(iEvent);
331       }
332       hlttree->Fill();
333       continue;
334     }
335
336     AliInfo(Form("processing event %d", iEvent));
337     fRunLoader->GetEvent(iEvent);
338
339     char fileName[256];
340     sprintf(fileName, "ESD_%d.%d_final.root", 
341             fRunLoader->GetHeader()->GetRun(), 
342             fRunLoader->GetHeader()->GetEventNrInRun());
343     if (!gSystem->AccessPathName(fileName)) continue;
344
345     // local reconstruction
346     if (!fRunLocalReconstruction.IsNull()) {
347       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
348         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
349       }
350     }
351
352     esd = new AliESD; hltesd = new AliESD;
353     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
354     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
355     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
356     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
357     if (gAlice) {
358       esd->SetMagneticField(gAlice->Field()->SolenoidField());
359       hltesd->SetMagneticField(gAlice->Field()->SolenoidField());
360     } else {
361       // ???
362     }
363
364     // vertex finder
365     if (fRunVertexFinder) {
366       if (!ReadESD(esd, "vertex")) {
367         if (!RunVertexFinder(esd)) {
368           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
369         }
370         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
371       }
372     }
373
374     // HLT tracking
375     if (!fRunTracking.IsNull()) {
376       if (fRunHLTTracking) {
377         hltesd->SetVertex(esd->GetVertex());
378         if (!RunHLTTracking(hltesd)) {
379           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
380         }
381       }
382     }
383
384     // barrel tracking
385     if (!fRunTracking.IsNull()) {
386       if (!ReadESD(esd, "tracking")) {
387         if (!RunTracking(esd)) {
388           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
389         }
390         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
391       }
392     }
393
394     // fill ESD
395     if (!fFillESD.IsNull()) {
396       if (!FillESD(esd, fFillESD)) {
397         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
398       }
399     }
400
401     // combined PID
402     AliESDpid::MakePID(esd);
403     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
404
405     // write ESD
406     tree->Fill();
407     // write HLT ESD
408     hlttree->Fill();
409
410     if (fCheckPointLevel > 0) WriteESD(esd, "final");
411     delete esd; delete hltesd;
412     esd = NULL; hltesd = NULL;
413   }
414
415   file->cd();
416   tree->Write();
417   hlttree->Write();
418   CleanUp(file, fileOld);
419
420   return kTRUE;
421 }
422
423
424 //_____________________________________________________________________________
425 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
426 {
427 // run the local reconstruction
428
429   TStopwatch stopwatch;
430   stopwatch.Start();
431
432   TString detStr = detectors;
433   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
434     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
435     AliReconstructor* reconstructor = GetReconstructor(iDet);
436     if (!reconstructor) continue;
437     if (reconstructor->HasLocalReconstruction()) continue;
438
439     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
440     TStopwatch stopwatchDet;
441     stopwatchDet.Start();
442     if (fRawReader) {
443       fRawReader->RewindEvents();
444       reconstructor->Reconstruct(fRunLoader, fRawReader);
445     } else {
446       reconstructor->Reconstruct(fRunLoader);
447     }
448     AliInfo(Form("execution time for %s:", fgkDetectorName[iDet]));
449     ToAliInfo(stopwatchDet.Print());
450   }
451
452   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
453     AliError(Form("the following detectors were not found: %s",
454                   detStr.Data()));
455     if (fStopOnError) return kFALSE;
456   }
457
458   AliInfo("execution time:");
459   ToAliInfo(stopwatch.Print());
460
461   return kTRUE;
462 }
463
464 //_____________________________________________________________________________
465 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
466 {
467 // run the local reconstruction
468
469   TStopwatch stopwatch;
470   stopwatch.Start();
471
472   TString detStr = detectors;
473   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
474     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
475     AliReconstructor* reconstructor = GetReconstructor(iDet);
476     if (!reconstructor) continue;
477     AliLoader* loader = fLoader[iDet];
478
479     // conversion of digits
480     if (fRawReader && reconstructor->HasDigitConversion()) {
481       AliInfo(Form("converting raw data digits into root objects for %s", 
482                    fgkDetectorName[iDet]));
483       TStopwatch stopwatchDet;
484       stopwatchDet.Start();
485       loader->LoadDigits("update");
486       loader->CleanDigits();
487       loader->MakeDigitsContainer();
488       TTree* digitsTree = loader->TreeD();
489       reconstructor->ConvertDigits(fRawReader, digitsTree);
490       loader->WriteDigits("OVERWRITE");
491       loader->UnloadDigits();
492       AliDebug(1, Form("execution time for %s:", fgkDetectorName[iDet]));
493       ToAliDebug(1, stopwatchDet.Print());
494     }
495
496     // local reconstruction
497     if (!reconstructor->HasLocalReconstruction()) continue;
498     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
499     TStopwatch stopwatchDet;
500     stopwatchDet.Start();
501     loader->LoadRecPoints("update");
502     loader->CleanRecPoints();
503     loader->MakeRecPointsContainer();
504     TTree* clustersTree = loader->TreeR();
505     if (fRawReader && !reconstructor->HasDigitConversion()) {
506       reconstructor->Reconstruct(fRawReader, clustersTree);
507     } else {
508       loader->LoadDigits("read");
509       TTree* digitsTree = loader->TreeD();
510       if (!digitsTree) {
511         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
512         if (fStopOnError) return kFALSE;
513       } else {
514         reconstructor->Reconstruct(digitsTree, clustersTree);
515       }
516       loader->UnloadDigits();
517     }
518     loader->WriteRecPoints("OVERWRITE");
519     loader->UnloadRecPoints();
520     AliDebug(1, Form("execution time for %s:", fgkDetectorName[iDet]));
521     ToAliDebug(1, stopwatchDet.Print());
522   }
523
524   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
525     AliError(Form("the following detectors were not found: %s",
526                   detStr.Data()));
527     if (fStopOnError) return kFALSE;
528   }
529
530   AliInfo("execution time:");
531   ToAliInfo(stopwatch.Print());
532
533   return kTRUE;
534 }
535
536 //_____________________________________________________________________________
537 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
538 {
539 // run the barrel tracking
540
541   TStopwatch stopwatch;
542   stopwatch.Start();
543
544   AliESDVertex* vertex = NULL;
545   Double_t vtxPos[3] = {0, 0, 0};
546   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
547   TArrayF mcVertex(3); 
548   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
549     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
550     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
551   }
552
553   if (fVertexer) {
554     AliInfo("running the ITS vertex finder");
555     if (fLoader[0]) fLoader[0]->LoadRecPoints();
556     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
557     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
558     if(!vertex){
559       AliWarning("Vertex not found");
560       vertex = new AliESDVertex();
561     }
562     else {
563       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
564     }
565
566   } else {
567     AliInfo("getting the primary vertex from MC");
568     vertex = new AliESDVertex(vtxPos, vtxErr);
569   }
570
571   if (vertex) {
572     vertex->GetXYZ(vtxPos);
573     vertex->GetSigmaXYZ(vtxErr);
574   } else {
575     AliWarning("no vertex reconstructed");
576     vertex = new AliESDVertex(vtxPos, vtxErr);
577   }
578   esd->SetVertex(vertex);
579   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
580     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
581   }  
582   delete vertex;
583
584   AliInfo("execution time:");
585   ToAliInfo(stopwatch.Print());
586
587   return kTRUE;
588 }
589
590 //_____________________________________________________________________________
591 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
592 {
593 // run the HLT barrel tracking
594
595   TStopwatch stopwatch;
596   stopwatch.Start();
597
598   if (!fRunLoader) {
599     AliError("Missing runLoader!");
600     return kFALSE;
601   }
602
603   AliInfo("running HLT tracking");
604
605   // Get a pointer to the HLT reconstructor
606   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
607   if (!reconstructor) return kFALSE;
608
609   // TPC + ITS
610   for (Int_t iDet = 1; iDet >= 0; iDet--) {
611     TString detName = fgkDetectorName[iDet];
612     AliDebug(1, Form("%s HLT tracking", detName.Data()));
613     reconstructor->SetOption(detName.Data());
614     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
615     if (!tracker) {
616       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
617       if (fStopOnError) return kFALSE;
618     }
619     Double_t vtxPos[3];
620     Double_t vtxErr[3]={0.005,0.005,0.010};
621     const AliESDVertex *vertex = esd->GetVertex();
622     vertex->GetXYZ(vtxPos);
623     tracker->SetVertex(vtxPos,vtxErr);
624     if(iDet != 1) {
625       fLoader[iDet]->LoadRecPoints("read");
626       TTree* tree = fLoader[iDet]->TreeR();
627       if (!tree) {
628         AliError(Form("Can't get the %s cluster tree", detName.Data()));
629         return kFALSE;
630       }
631       tracker->LoadClusters(tree);
632     }
633     if (tracker->Clusters2Tracks(esd) != 0) {
634       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
635       return kFALSE;
636     }
637     if(iDet != 1) {
638       tracker->UnloadClusters();
639     }
640     delete tracker;
641   }
642
643   AliInfo("execution time:");
644   ToAliInfo(stopwatch.Print());
645
646   return kTRUE;
647 }
648
649 //_____________________________________________________________________________
650 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
651 {
652 // run the barrel tracking
653
654   TStopwatch stopwatch;
655   stopwatch.Start();
656
657   AliInfo("running tracking");
658
659   // pass 1: TPC + ITS inwards
660   for (Int_t iDet = 1; iDet >= 0; iDet--) {
661     if (!fTracker[iDet]) continue;
662     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
663
664     // load clusters
665     fLoader[iDet]->LoadRecPoints("read");
666     TTree* tree = fLoader[iDet]->TreeR();
667     if (!tree) {
668       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
669       return kFALSE;
670     }
671     fTracker[iDet]->LoadClusters(tree);
672
673     // run tracking
674     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
675       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
676       return kFALSE;
677     }
678     if (fCheckPointLevel > 1) {
679       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
680     }
681     // preliminary PID in TPC needed by the ITS tracker
682     if (iDet == 1) {
683       GetReconstructor(1)->FillESD(fRunLoader, esd);
684       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
685       AliESDpid::MakePID(esd);
686     }
687   }
688
689   // pass 2: ALL backwards
690   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
691     if (!fTracker[iDet]) continue;
692     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
693
694     // load clusters
695     if (iDet > 1) {     // all except ITS, TPC
696       TTree* tree = NULL;
697       if (iDet == 3) {   // TOF
698         fLoader[iDet]->LoadDigits("read");
699         tree = fLoader[iDet]->TreeD();
700       } else {
701         fLoader[iDet]->LoadRecPoints("read");
702         tree = fLoader[iDet]->TreeR();
703       }
704       if (!tree) {
705         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
706         return kFALSE;
707       }
708       fTracker[iDet]->LoadClusters(tree);
709     }
710
711     // run tracking
712     if (fTracker[iDet]->PropagateBack(esd) != 0) {
713       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
714       return kFALSE;
715     }
716     if (fCheckPointLevel > 1) {
717       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
718     }
719
720     // unload clusters
721     if (iDet > 2) {     // all except ITS, TPC, TRD
722       fTracker[iDet]->UnloadClusters();
723       if (iDet == 3) {   // TOF
724         fLoader[iDet]->UnloadDigits();
725       } else {
726         fLoader[iDet]->UnloadRecPoints();
727       }
728     }
729   }
730
731   // pass 3: TRD + TPC + ITS refit inwards
732   for (Int_t iDet = 2; iDet >= 0; iDet--) {
733     if (!fTracker[iDet]) continue;
734     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
735
736     // run tracking
737     if (fTracker[iDet]->RefitInward(esd) != 0) {
738       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
739       return kFALSE;
740     }
741     if (fCheckPointLevel > 1) {
742       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
743     }
744
745     // unload clusters
746     fTracker[iDet]->UnloadClusters();
747     fLoader[iDet]->UnloadRecPoints();
748   }
749
750   AliInfo("execution time:");
751   ToAliInfo(stopwatch.Print());
752
753   return kTRUE;
754 }
755
756 //_____________________________________________________________________________
757 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
758 {
759 // fill the event summary data
760
761   TStopwatch stopwatch;
762   stopwatch.Start();
763   AliInfo("filling ESD");
764
765   TString detStr = detectors;
766   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
767     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
768     AliReconstructor* reconstructor = GetReconstructor(iDet);
769     if (!reconstructor) continue;
770
771     if (!ReadESD(esd, fgkDetectorName[iDet])) {
772       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
773       TTree* clustersTree = NULL;
774       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
775         fLoader[iDet]->LoadRecPoints("read");
776         clustersTree = fLoader[iDet]->TreeR();
777         if (!clustersTree) {
778           AliError(Form("Can't get the %s clusters tree", 
779                         fgkDetectorName[iDet]));
780           if (fStopOnError) return kFALSE;
781         }
782       }
783       if (fRawReader && !reconstructor->HasDigitConversion()) {
784         reconstructor->FillESD(fRawReader, clustersTree, esd);
785       } else {
786         TTree* digitsTree = NULL;
787         if (fLoader[iDet]) {
788           fLoader[iDet]->LoadDigits("read");
789           digitsTree = fLoader[iDet]->TreeD();
790           if (!digitsTree) {
791             AliError(Form("Can't get the %s digits tree", 
792                           fgkDetectorName[iDet]));
793             if (fStopOnError) return kFALSE;
794           }
795         }
796         reconstructor->FillESD(digitsTree, clustersTree, esd);
797         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
798       }
799       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
800         fLoader[iDet]->UnloadRecPoints();
801       }
802
803       if (fRawReader) {
804         reconstructor->FillESD(fRunLoader, fRawReader, esd);
805       } else {
806         reconstructor->FillESD(fRunLoader, esd);
807       }
808       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
809     }
810   }
811
812   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
813     AliError(Form("the following detectors were not found: %s", 
814                   detStr.Data()));
815     if (fStopOnError) return kFALSE;
816   }
817
818   AliInfo("execution time:");
819   ToAliInfo(stopwatch.Print());
820
821   return kTRUE;
822 }
823
824
825 //_____________________________________________________________________________
826 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
827 {
828 // check whether detName is contained in detectors
829 // if yes, it is removed from detectors
830
831   // check if all detectors are selected
832   if ((detectors.CompareTo("ALL") == 0) ||
833       detectors.BeginsWith("ALL ") ||
834       detectors.EndsWith(" ALL") ||
835       detectors.Contains(" ALL ")) {
836     detectors = "ALL";
837     return kTRUE;
838   }
839
840   // search for the given detector
841   Bool_t result = kFALSE;
842   if ((detectors.CompareTo(detName) == 0) ||
843       detectors.BeginsWith(detName+" ") ||
844       detectors.EndsWith(" "+detName) ||
845       detectors.Contains(" "+detName+" ")) {
846     detectors.ReplaceAll(detName, "");
847     result = kTRUE;
848   }
849
850   // clean up the detectors string
851   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
852   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
853   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
854
855   return result;
856 }
857
858 //_____________________________________________________________________________
859 Bool_t AliReconstruction::InitRunLoader()
860 {
861 // get or create the run loader
862
863   if (gAlice) delete gAlice;
864   gAlice = NULL;
865
866   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
867     // load all base libraries to get the loader classes
868     TString libs = gSystem->GetLibraries();
869     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
870       TString detName = fgkDetectorName[iDet];
871       if (detName == "HLT") continue;
872       if (libs.Contains("lib" + detName + "base.so")) continue;
873       gSystem->Load("lib" + detName + "base.so");
874     }
875     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
876     if (!fRunLoader) {
877       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
878       CleanUp();
879       return kFALSE;
880     }
881     fRunLoader->CdGAFile();
882     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
883       if (fRunLoader->LoadgAlice() == 0) {
884         gAlice = fRunLoader->GetAliRun();
885         AliTracker::SetFieldMap(gAlice->Field());
886       }
887     }
888     if (!gAlice && !fRawReader) {
889       AliError(Form("no gAlice object found in file %s",
890                     fGAliceFileName.Data()));
891       CleanUp();
892       return kFALSE;
893     }
894
895   } else {               // galice.root does not exist
896     if (!fRawReader) {
897       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
898       CleanUp();
899       return kFALSE;
900     }
901     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
902                                     AliConfig::GetDefaultEventFolderName(),
903                                     "recreate");
904     if (!fRunLoader) {
905       AliError(Form("could not create run loader in file %s", 
906                     fGAliceFileName.Data()));
907       CleanUp();
908       return kFALSE;
909     }
910     fRunLoader->MakeTree("E");
911     Int_t iEvent = 0;
912     while (fRawReader->NextEvent()) {
913       fRunLoader->SetEventNumber(iEvent);
914       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
915                                      iEvent, iEvent);
916       fRunLoader->MakeTree("H");
917       fRunLoader->TreeE()->Fill();
918       iEvent++;
919     }
920     fRawReader->RewindEvents();
921     fRunLoader->WriteHeader("OVERWRITE");
922     fRunLoader->CdGAFile();
923     fRunLoader->Write(0, TObject::kOverwrite);
924 //    AliTracker::SetFieldMap(???);
925   }
926
927   return kTRUE;
928 }
929
930 //_____________________________________________________________________________
931 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
932 {
933 // get the reconstructor object and the loader for a detector
934
935   if (fReconstructor[iDet]) return fReconstructor[iDet];
936
937   // load the reconstructor object
938   TPluginManager* pluginManager = gROOT->GetPluginManager();
939   TString detName = fgkDetectorName[iDet];
940   TString recName = "Ali" + detName + "Reconstructor";
941   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
942
943   if (detName == "HLT") {
944     if (!gROOT->GetClass("AliLevel3")) {
945       gSystem->Load("libAliL3Src.so");
946       gSystem->Load("libAliL3Misc.so");
947       gSystem->Load("libAliL3Hough.so");
948       gSystem->Load("libAliL3Comp.so");
949     }
950   }
951
952   AliReconstructor* reconstructor = NULL;
953   // first check if a plugin is defined for the reconstructor
954   TPluginHandler* pluginHandler = 
955     pluginManager->FindHandler("AliReconstructor", detName);
956   // if not, add a plugin for it
957   if (!pluginHandler) {
958     AliDebug(1, Form("defining plugin for %s", recName.Data()));
959     TString libs = gSystem->GetLibraries();
960     if (libs.Contains("lib" + detName + "base.so") ||
961         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
962       pluginManager->AddHandler("AliReconstructor", detName, 
963                                 recName, detName + "rec", recName + "()");
964     } else {
965       pluginManager->AddHandler("AliReconstructor", detName, 
966                                 recName, detName, recName + "()");
967     }
968     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
969   }
970   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
971     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
972   }
973   if (reconstructor) {
974     TObject* obj = fOptions.FindObject(detName.Data());
975     if (obj) reconstructor->SetOption(obj->GetTitle());
976     reconstructor->Init(fRunLoader);
977     fReconstructor[iDet] = reconstructor;
978   }
979
980   // get or create the loader
981   if (detName != "HLT") {
982     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
983     if (!fLoader[iDet]) {
984       AliConfig::Instance()
985         ->CreateDetectorFolders(fRunLoader->GetEventFolder(), 
986                                 detName, detName);
987       // first check if a plugin is defined for the loader
988       TPluginHandler* pluginHandler = 
989         pluginManager->FindHandler("AliLoader", detName);
990       // if not, add a plugin for it
991       if (!pluginHandler) {
992         TString loaderName = "Ali" + detName + "Loader";
993         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
994         pluginManager->AddHandler("AliLoader", detName, 
995                                   loaderName, detName + "base", 
996                                   loaderName + "(const char*, TFolder*)");
997         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
998       }
999       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1000         fLoader[iDet] = 
1001           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), 
1002                                                  fRunLoader->GetEventFolder());
1003       }
1004       if (!fLoader[iDet]) {   // use default loader
1005         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1006       }
1007       if (!fLoader[iDet]) {
1008         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1009         if (fStopOnError) return NULL;
1010       } else {
1011         fRunLoader->AddLoader(fLoader[iDet]);
1012         fRunLoader->CdGAFile();
1013         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1014         fRunLoader->Write(0, TObject::kOverwrite);
1015       }
1016     }
1017   }
1018       
1019   return reconstructor;
1020 }
1021
1022 //_____________________________________________________________________________
1023 Bool_t AliReconstruction::CreateVertexer()
1024 {
1025 // create the vertexer
1026
1027   fVertexer = NULL;
1028   AliReconstructor* itsReconstructor = GetReconstructor(0);
1029   if (itsReconstructor) {
1030     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1031   }
1032   if (!fVertexer) {
1033     AliWarning("couldn't create a vertexer for ITS");
1034     if (fStopOnError) return kFALSE;
1035   }
1036
1037   return kTRUE;
1038 }
1039
1040 //_____________________________________________________________________________
1041 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1042 {
1043 // create the trackers
1044
1045   TString detStr = detectors;
1046   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1047     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1048     AliReconstructor* reconstructor = GetReconstructor(iDet);
1049     if (!reconstructor) continue;
1050     TString detName = fgkDetectorName[iDet];
1051     if (detName == "HLT") {
1052       fRunHLTTracking = kTRUE;
1053       continue;
1054     }
1055
1056     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1057     if (!fTracker[iDet] && (iDet < 7)) {
1058       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1059       if (fStopOnError) return kFALSE;
1060     }
1061   }
1062
1063   return kTRUE;
1064 }
1065
1066 //_____________________________________________________________________________
1067 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1068 {
1069 // delete trackers and the run loader and close and delete the file
1070
1071   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1072     delete fReconstructor[iDet];
1073     fReconstructor[iDet] = NULL;
1074     fLoader[iDet] = NULL;
1075     delete fTracker[iDet];
1076     fTracker[iDet] = NULL;
1077   }
1078   delete fVertexer;
1079   fVertexer = NULL;
1080
1081   delete fRunLoader;
1082   fRunLoader = NULL;
1083   delete fRawReader;
1084   fRawReader = NULL;
1085
1086   if (file) {
1087     file->Close();
1088     delete file;
1089   }
1090
1091   if (fileOld) {
1092     fileOld->Close();
1093     delete fileOld;
1094     gSystem->Unlink("AliESDs.old.root");
1095   }
1096 }
1097
1098
1099 //_____________________________________________________________________________
1100 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1101 {
1102 // read the ESD event from a file
1103
1104   if (!esd) return kFALSE;
1105   char fileName[256];
1106   sprintf(fileName, "ESD_%d.%d_%s.root", 
1107           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1108   if (gSystem->AccessPathName(fileName)) return kFALSE;
1109
1110   AliDebug(1, Form("reading ESD from file %s", fileName));
1111   TFile* file = TFile::Open(fileName);
1112   if (!file || !file->IsOpen()) {
1113     AliError(Form("opening %s failed", fileName));
1114     delete file;
1115     return kFALSE;
1116   }
1117
1118   gROOT->cd();
1119   delete esd;
1120   esd = (AliESD*) file->Get("ESD");
1121   file->Close();
1122   delete file;
1123   return kTRUE;
1124 }
1125
1126 //_____________________________________________________________________________
1127 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1128 {
1129 // write the ESD event to a file
1130
1131   if (!esd) return;
1132   char fileName[256];
1133   sprintf(fileName, "ESD_%d.%d_%s.root", 
1134           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1135
1136   AliDebug(1, Form("writing ESD to file %s", fileName));
1137   TFile* file = TFile::Open(fileName, "recreate");
1138   if (!file || !file->IsOpen()) {
1139     AliError(Form("opening %s failed", fileName));
1140   } else {
1141     esd->Write("ESD");
1142     file->Close();
1143   }
1144   delete file;
1145 }