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