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