417388ef6409315e938868eb23038e84f29fc0f6
[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.SetUniformFieldTracking(kFALSE); ) //
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 #include <TGeoManager.h>
115 #include <TLorentzVector.h>
116
117 #include "AliReconstruction.h"
118 #include "AliReconstructor.h"
119 #include "AliLog.h"
120 #include "AliRunLoader.h"
121 #include "AliRun.h"
122 #include "AliRawReaderFile.h"
123 #include "AliRawReaderDate.h"
124 #include "AliRawReaderRoot.h"
125 #include "AliESD.h"
126 #include "AliESDfriend.h"
127 #include "AliESDVertex.h"
128 #include "AliTracker.h"
129 #include "AliVertexer.h"
130 #include "AliVertexerTracks.h"
131 #include "AliHeader.h"
132 #include "AliGenEventHeader.h"
133 #include "AliPID.h"
134 #include "AliESDpid.h"
135 #include "AliESDtrack.h"
136
137 #include "AliRunTag.h"
138 #include "AliDetectorTag.h"
139 #include "AliEventTag.h"
140
141 #include "AliTrackPointArray.h"
142 #include "AliCDBManager.h"
143 #include "AliCDBEntry.h"
144 #include "AliAlignObj.h"
145
146 #include "AliCentralTrigger.h"
147 #include "AliCTPRawStream.h"
148
149 ClassImp(AliReconstruction)
150
151
152 //_____________________________________________________________________________
153 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
154
155 //_____________________________________________________________________________
156 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
157                                      const char* name, const char* title) :
158   TNamed(name, title),
159
160   fUniformField(kTRUE),
161   fRunVertexFinder(kTRUE),
162   fRunHLTTracking(kFALSE),
163   fStopOnError(kFALSE),
164   fWriteAlignmentData(kFALSE),
165   fWriteESDfriend(kFALSE),
166   fFillTriggerESD(kTRUE),
167
168   fRunLocalReconstruction("ALL"),
169   fRunTracking("ALL"),
170   fFillESD("ALL"),
171   fGAliceFileName(gAliceFilename),
172   fInput(""),
173   fFirstEvent(0),
174   fLastEvent(-1),
175   fCheckPointLevel(0),
176   fOptions(),
177   fLoadAlignFromCDB(kTRUE),
178   fLoadAlignData("ALL"),
179
180   fRunLoader(NULL),
181   fRawReader(NULL),
182
183   fVertexer(NULL),
184
185   fAlignObjArray(NULL),
186   fCDBUri(cdbUri)
187 {
188 // create reconstruction object with default parameters
189   
190   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
191     fReconstructor[iDet] = NULL;
192     fLoader[iDet] = NULL;
193     fTracker[iDet] = NULL;
194   }
195   AliPID pid;
196 }
197
198 //_____________________________________________________________________________
199 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
200   TNamed(rec),
201
202   fUniformField(rec.fUniformField),
203   fRunVertexFinder(rec.fRunVertexFinder),
204   fRunHLTTracking(rec.fRunHLTTracking),
205   fStopOnError(rec.fStopOnError),
206   fWriteAlignmentData(rec.fWriteAlignmentData),
207   fWriteESDfriend(rec.fWriteESDfriend),
208   fFillTriggerESD(rec.fFillTriggerESD),
209
210   fRunLocalReconstruction(rec.fRunLocalReconstruction),
211   fRunTracking(rec.fRunTracking),
212   fFillESD(rec.fFillESD),
213   fGAliceFileName(rec.fGAliceFileName),
214   fInput(rec.fInput),
215   fFirstEvent(rec.fFirstEvent),
216   fLastEvent(rec.fLastEvent),
217   fCheckPointLevel(0),
218   fOptions(),
219   fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
220   fLoadAlignData(rec.fLoadAlignData),
221
222   fRunLoader(NULL),
223   fRawReader(NULL),
224
225   fVertexer(NULL),
226
227   fAlignObjArray(rec.fAlignObjArray),
228   fCDBUri(rec.fCDBUri)
229 {
230 // copy constructor
231
232   for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
233     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
234   }
235   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
236     fReconstructor[iDet] = NULL;
237     fLoader[iDet] = NULL;
238     fTracker[iDet] = NULL;
239   }
240 }
241
242 //_____________________________________________________________________________
243 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
244 {
245 // assignment operator
246
247   this->~AliReconstruction();
248   new(this) AliReconstruction(rec);
249   return *this;
250 }
251
252 //_____________________________________________________________________________
253 AliReconstruction::~AliReconstruction()
254 {
255 // clean up
256
257   CleanUp();
258   fOptions.Delete();
259 }
260
261 //_____________________________________________________________________________
262 void AliReconstruction::InitCDBStorage()
263 {
264 // activate a default CDB storage
265 // First check if we have any CDB storage set, because it is used 
266 // to retrieve the calibration and alignment constants
267
268   AliCDBManager* man = AliCDBManager::Instance();
269   if (!man->IsDefaultStorageSet())
270   {
271     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
272     AliWarning("Default CDB storage not yet set");
273     AliWarning(Form("Using default storage declared in AliSimulation: %s",fCDBUri.Data()));
274     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
275     SetDefaultStorage(fCDBUri);
276     
277     Int_t cdbRun = AliCDBManager::Instance()->GetRun();
278     if(cdbRun == -1){
279         AliWarning("AliCDBManager's run number temporarily set to 0!!");
280         AliCDBManager::Instance()->SetRun(0);
281     }
282   }  
283   
284 }
285
286 //_____________________________________________________________________________
287 void AliReconstruction::SetDefaultStorage(const char* uri) {
288 // activate a default CDB storage 
289
290    AliCDBManager::Instance()->SetDefaultStorage(uri);
291
292 }
293
294 //_____________________________________________________________________________
295 void AliReconstruction::SetSpecificStorage(const char* detName, const char* uri) {
296 // activate a detector-specific CDB storage 
297
298    AliCDBManager::Instance()->SetSpecificStorage(detName, uri);
299
300 }
301
302 //_____________________________________________________________________________
303 Bool_t AliReconstruction::SetRunNumber()
304 {
305   // The method is called in Run() in order
306   // to set a correct run number.
307   // In case of raw data reconstruction the
308   // run number is taken from the raw data header
309
310   if(AliCDBManager::Instance()->GetRun() < 0) {
311     if (!fRunLoader) {
312       AliError("No run loader is found !"); 
313       return kFALSE;
314     }
315     // read run number from gAlice
316     AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
317     AliInfo(Form("Run number: %d",AliCDBManager::Instance()->GetRun()));
318   }
319   return kTRUE;
320 }
321
322 //_____________________________________________________________________________
323 Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
324 {
325   // Read collection of alignment objects (AliAlignObj derived) saved
326   // in the TClonesArray ClArrayName and apply them to the geometry
327   // manager singleton.
328   //
329   alObjArray->Sort();
330   Int_t nvols = alObjArray->GetEntriesFast();
331
332   Bool_t flag = kTRUE;
333
334   for(Int_t j=0; j<nvols; j++)
335     {
336       AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
337       if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
338     }
339
340   if (AliDebugLevelClass() >= 1) {
341     gGeoManager->CheckOverlaps(20);
342     TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
343     if(ovexlist->GetEntriesFast()){  
344       AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
345    }
346   }
347
348   return flag;
349
350 }
351
352 //_____________________________________________________________________________
353 Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
354 {
355   // Fills array of single detector's alignable objects from CDB
356   
357   AliDebug(2, Form("Loading alignment data for detector: %s",detName));
358   
359   AliCDBEntry *entry;
360         
361   AliCDBPath path(detName,"Align","Data");
362         
363   entry=AliCDBManager::Instance()->Get(path.GetPath());
364   if(!entry){ 
365         AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
366         return kFALSE;
367   }
368   entry->SetOwner(1);
369   TClonesArray *alignArray = (TClonesArray*) entry->GetObject();        
370   alignArray->SetOwner(0);
371   AliDebug(2,Form("Found %d alignment objects for %s",
372                         alignArray->GetEntries(),detName));
373
374   AliAlignObj *alignObj=0;
375   TIter iter(alignArray);
376         
377   // loop over align objects in detector
378   while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
379         fAlignObjArray->Add(alignObj);
380   }
381   // delete entry --- Don't delete, it is cached!
382         
383   AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
384   return kTRUE;
385
386 }
387
388 //_____________________________________________________________________________
389 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
390 {
391   // Read the alignment objects from CDB.
392   // Each detector is supposed to have the
393   // alignment objects in DET/Align/Data CDB path.
394   // All the detector objects are then collected,
395   // sorted by geometry level (starting from ALIC) and
396   // then applied to the TGeo geometry.
397   // Finally an overlaps check is performed.
398
399   // Load alignment data from CDB and fill fAlignObjArray 
400   if(fLoadAlignFromCDB){
401         if(!fAlignObjArray) fAlignObjArray = new TObjArray();
402         
403         //fAlignObjArray->RemoveAll(); 
404         fAlignObjArray->Clear();        
405         fAlignObjArray->SetOwner(0);
406  
407         TString detStr = detectors;
408         TString dataNotLoaded="";
409         TString dataLoaded="";
410   
411         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
412           if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
413           if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
414             dataNotLoaded += fgkDetectorName[iDet];
415             dataNotLoaded += " ";
416           } else {
417             dataLoaded += fgkDetectorName[iDet];
418             dataLoaded += " ";
419           }
420         } // end loop over detectors
421   
422         if ((detStr.CompareTo("ALL") == 0)) detStr = "";
423         dataNotLoaded += detStr;
424         AliInfo(Form("Alignment data loaded for: %s",
425                           dataLoaded.Data()));
426         AliInfo(Form("Didn't/couldn't load alignment data for: %s",
427                           dataNotLoaded.Data()));
428   } // fLoadAlignFromCDB flag
429  
430   // Check if the array with alignment objects was
431   // provided by the user. If yes, apply the objects
432   // to the present TGeo geometry
433   if (fAlignObjArray) {
434     if (gGeoManager && gGeoManager->IsClosed()) {
435       if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
436         AliError("The misalignment of one or more volumes failed!"
437                  "Compare the list of simulated detectors and the list of detector alignment data!");
438         return kFALSE;
439       }
440     }
441     else {
442       AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
443       return kFALSE;
444     }
445   }
446
447   return kTRUE;
448 }
449
450 //_____________________________________________________________________________
451 void AliReconstruction::SetGAliceFile(const char* fileName)
452 {
453 // set the name of the galice file
454
455   fGAliceFileName = fileName;
456 }
457
458 //_____________________________________________________________________________
459 void AliReconstruction::SetOption(const char* detector, const char* option)
460 {
461 // set options for the reconstruction of a detector
462
463   TObject* obj = fOptions.FindObject(detector);
464   if (obj) fOptions.Remove(obj);
465   fOptions.Add(new TNamed(detector, option));
466 }
467
468
469 //_____________________________________________________________________________
470 Bool_t AliReconstruction::Run(const char* input,
471                               Int_t firstEvent, Int_t lastEvent)
472 {
473 // run the reconstruction
474
475   InitCDBStorage();
476
477   // set the input
478   if (!input) input = fInput.Data();
479   TString fileName(input);
480   if (fileName.EndsWith("/")) {
481     fRawReader = new AliRawReaderFile(fileName);
482   } else if (fileName.EndsWith(".root")) {
483     fRawReader = new AliRawReaderRoot(fileName);
484   } else if (!fileName.IsNull()) {
485     fRawReader = new AliRawReaderDate(fileName);
486     fRawReader->SelectEvents(7);
487   }
488
489   // get the run loader
490   if (!InitRunLoader()) return kFALSE;
491
492   // Set run number in CDBManager (if it is not already set by the user)
493   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
494
495   // Import ideal TGeo geometry and apply misalignment
496   if (!gGeoManager) {
497     TString geom(gSystem->DirName(fGAliceFileName));
498     geom += "/geometry.root";
499     TGeoManager::Import(geom.Data());
500     if (!gGeoManager) if (fStopOnError) return kFALSE;
501   }
502   if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
503
504   // Temporary fix by A.Gheata
505   // Could be removed with the next Root version (>5.11)
506   if (gGeoManager) {
507     TIter next(gGeoManager->GetListOfVolumes());
508     TGeoVolume *vol;
509     while ((vol = (TGeoVolume *)next())) {
510       if (vol->GetVoxels()) {
511         if (vol->GetVoxels()->NeedRebuild()) {
512           vol->GetVoxels()->Voxelize();
513           vol->FindOverlaps();
514         }
515       }
516     }
517   }
518
519   // local reconstruction
520   if (!fRunLocalReconstruction.IsNull()) {
521     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
522       if (fStopOnError) {CleanUp(); return kFALSE;}
523     }
524   }
525 //  if (!fRunVertexFinder && fRunTracking.IsNull() && 
526 //      fFillESD.IsNull()) return kTRUE;
527
528   // get vertexer
529   if (fRunVertexFinder && !CreateVertexer()) {
530     if (fStopOnError) {
531       CleanUp(); 
532       return kFALSE;
533     }
534   }
535
536   // get trackers
537   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
538     if (fStopOnError) {
539       CleanUp(); 
540       return kFALSE;
541     }      
542   }
543
544
545   TStopwatch stopwatch;
546   stopwatch.Start();
547
548   // get the possibly already existing ESD file and tree
549   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
550   TFile* fileOld = NULL;
551   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
552   if (!gSystem->AccessPathName("AliESDs.root")){
553     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
554     fileOld = TFile::Open("AliESDs.old.root");
555     if (fileOld && fileOld->IsOpen()) {
556       treeOld = (TTree*) fileOld->Get("esdTree");
557       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
558       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
559       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
560     }
561   }
562
563   // create the ESD output file and tree
564   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
565   if (!file->IsOpen()) {
566     AliError("opening AliESDs.root failed");
567     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}    
568   }
569   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
570   tree->Branch("ESD", "AliESD", &esd);
571   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
572   hlttree->Branch("ESD", "AliESD", &hltesd);
573   delete esd; delete hltesd;
574   esd = NULL; hltesd = NULL;
575
576   // create the file and tree with ESD additions
577   TFile *filef=0; TTree *treef=0; AliESDfriend *esdf=0;
578   if (fWriteESDfriend) {
579      filef = TFile::Open("AliESDfriends.root", "RECREATE");
580      if (!filef->IsOpen()) {
581         AliError("opening AliESDfriends.root failed");
582      }
583      treef = new TTree("esdFriendTree", "Tree with ESD friends");
584      treef->Branch("ESDfriend", "AliESDfriend", &esdf);
585   }
586
587   gROOT->cd();
588
589   AliVertexerTracks tVertexer;
590
591   // loop over events
592   if (fRawReader) fRawReader->RewindEvents();
593   
594   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
595     if (fRawReader) fRawReader->NextEvent();
596     if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
597       // copy old ESD to the new one
598       if (treeOld) {
599         treeOld->SetBranchAddress("ESD", &esd);
600         treeOld->GetEntry(iEvent);
601       }
602       tree->Fill();
603       if (hlttreeOld) {
604         hlttreeOld->SetBranchAddress("ESD", &hltesd);
605         hlttreeOld->GetEntry(iEvent);
606       }
607       hlttree->Fill();
608       continue;
609     }
610
611     AliInfo(Form("processing event %d", iEvent));
612     fRunLoader->GetEvent(iEvent);
613
614     char fileName[256];
615     sprintf(fileName, "ESD_%d.%d_final.root", 
616             fRunLoader->GetHeader()->GetRun(), 
617             fRunLoader->GetHeader()->GetEventNrInRun());
618     if (!gSystem->AccessPathName(fileName)) continue;
619
620     // local reconstruction
621     if (!fRunLocalReconstruction.IsNull()) {
622       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
623         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
624       }
625     }
626
627     esd = new AliESD; hltesd = new AliESD;
628     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
629     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
630     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
631     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
632     if (gAlice) {
633       esd->SetMagneticField(AliTracker::GetBz());
634       hltesd->SetMagneticField(AliTracker::GetBz());
635     } else {
636       // ???
637     }
638
639     // vertex finder
640     if (fRunVertexFinder) {
641       if (!ReadESD(esd, "vertex")) {
642         if (!RunVertexFinder(esd)) {
643           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
644         }
645         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
646       }
647     }
648
649     // HLT tracking
650     if (!fRunTracking.IsNull()) {
651       if (fRunHLTTracking) {
652         hltesd->SetVertex(esd->GetVertex());
653         if (!RunHLTTracking(hltesd)) {
654           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
655         }
656       }
657     }
658
659     // barrel tracking
660     if (!fRunTracking.IsNull()) {
661       if (!ReadESD(esd, "tracking")) {
662         if (!RunTracking(esd)) {
663           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
664         }
665         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
666       }
667     }
668
669     // fill ESD
670     if (!fFillESD.IsNull()) {
671       if (!FillESD(esd, fFillESD)) {
672         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
673       }
674     }
675
676     // combined PID
677     AliESDpid::MakePID(esd);
678     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
679
680     if (fFillTriggerESD) {
681       if (!ReadESD(esd, "trigger")) {
682         if (!FillTriggerESD(esd)) {
683           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
684         }
685         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
686       }
687     }
688
689     esd->SetPrimaryVertex(tVertexer.FindPrimaryVertex(esd));
690
691     // write ESD
692     tree->Fill();
693     // write HLT ESD
694     hlttree->Fill();
695
696     // write ESD friend
697     if (fWriteESDfriend) {
698        esdf=new AliESDfriend();
699        esd->GetESDfriend(esdf);
700        treef->Fill();  
701     }
702
703     if (fCheckPointLevel > 0)  WriteESD(esd, "final"); 
704  
705     delete esd; delete hltesd;
706     esd = NULL; hltesd = NULL;
707   }
708
709   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
710                stopwatch.RealTime(),stopwatch.CpuTime()));
711
712   file->cd();
713   tree->Write();
714   hlttree->Write();
715
716   if (fWriteESDfriend) {
717      filef->cd();
718      treef->Write(); delete treef; filef->Close(); delete filef; 
719   }
720
721   // Create tags for the events in the ESD tree (the ESD tree is always present)
722   // In case of empty events the tags will contain dummy values
723   CreateTag(file);
724   CleanUp(file, fileOld);
725
726   return kTRUE;
727 }
728
729
730 //_____________________________________________________________________________
731 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
732 {
733 // run the local reconstruction
734
735   TStopwatch stopwatch;
736   stopwatch.Start();
737
738   TString detStr = detectors;
739   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
740     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
741     AliReconstructor* reconstructor = GetReconstructor(iDet);
742     if (!reconstructor) continue;
743     if (reconstructor->HasLocalReconstruction()) continue;
744
745     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
746     TStopwatch stopwatchDet;
747     stopwatchDet.Start();
748     if (fRawReader) {
749       fRawReader->RewindEvents();
750       reconstructor->Reconstruct(fRunLoader, fRawReader);
751     } else {
752       reconstructor->Reconstruct(fRunLoader);
753     }
754     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
755                  fgkDetectorName[iDet],
756                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
757   }
758
759   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
760     AliError(Form("the following detectors were not found: %s",
761                   detStr.Data()));
762     if (fStopOnError) return kFALSE;
763   }
764
765   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
766                stopwatch.RealTime(),stopwatch.CpuTime()));
767
768   return kTRUE;
769 }
770
771 //_____________________________________________________________________________
772 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
773 {
774 // run the local reconstruction
775
776   TStopwatch stopwatch;
777   stopwatch.Start();
778
779   TString detStr = detectors;
780   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
781     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
782     AliReconstructor* reconstructor = GetReconstructor(iDet);
783     if (!reconstructor) continue;
784     AliLoader* loader = fLoader[iDet];
785
786     // conversion of digits
787     if (fRawReader && reconstructor->HasDigitConversion()) {
788       AliInfo(Form("converting raw data digits into root objects for %s", 
789                    fgkDetectorName[iDet]));
790       TStopwatch stopwatchDet;
791       stopwatchDet.Start();
792       loader->LoadDigits("update");
793       loader->CleanDigits();
794       loader->MakeDigitsContainer();
795       TTree* digitsTree = loader->TreeD();
796       reconstructor->ConvertDigits(fRawReader, digitsTree);
797       loader->WriteDigits("OVERWRITE");
798       loader->UnloadDigits();
799       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
800                    fgkDetectorName[iDet],
801                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
802     }
803
804     // local reconstruction
805     if (!reconstructor->HasLocalReconstruction()) continue;
806     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
807     TStopwatch stopwatchDet;
808     stopwatchDet.Start();
809     loader->LoadRecPoints("update");
810     loader->CleanRecPoints();
811     loader->MakeRecPointsContainer();
812     TTree* clustersTree = loader->TreeR();
813     if (fRawReader && !reconstructor->HasDigitConversion()) {
814       reconstructor->Reconstruct(fRawReader, clustersTree);
815     } else {
816       loader->LoadDigits("read");
817       TTree* digitsTree = loader->TreeD();
818       if (!digitsTree) {
819         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
820         if (fStopOnError) return kFALSE;
821       } else {
822         reconstructor->Reconstruct(digitsTree, clustersTree);
823       }
824       loader->UnloadDigits();
825     }
826     loader->WriteRecPoints("OVERWRITE");
827     loader->UnloadRecPoints();
828     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
829                     fgkDetectorName[iDet],
830                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
831   }
832
833   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
834     AliError(Form("the following detectors were not found: %s",
835                   detStr.Data()));
836     if (fStopOnError) return kFALSE;
837   }
838   
839   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
840                stopwatch.RealTime(),stopwatch.CpuTime()));
841
842   return kTRUE;
843 }
844
845 //_____________________________________________________________________________
846 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
847 {
848 // run the barrel tracking
849
850   TStopwatch stopwatch;
851   stopwatch.Start();
852
853   AliESDVertex* vertex = NULL;
854   Double_t vtxPos[3] = {0, 0, 0};
855   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
856   TArrayF mcVertex(3); 
857   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
858     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
859     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
860   }
861
862   if (fVertexer) {
863     AliInfo("running the ITS vertex finder");
864     if (fLoader[0]) fLoader[0]->LoadRecPoints();
865     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
866     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
867     if(!vertex){
868       AliWarning("Vertex not found");
869       vertex = new AliESDVertex();
870       vertex->SetName("default");
871     }
872     else {
873       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
874       vertex->SetName("reconstructed");
875     }
876
877   } else {
878     AliInfo("getting the primary vertex from MC");
879     vertex = new AliESDVertex(vtxPos, vtxErr);
880   }
881
882   if (vertex) {
883     vertex->GetXYZ(vtxPos);
884     vertex->GetSigmaXYZ(vtxErr);
885   } else {
886     AliWarning("no vertex reconstructed");
887     vertex = new AliESDVertex(vtxPos, vtxErr);
888   }
889   esd->SetVertex(vertex);
890   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
891     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
892   }  
893   delete vertex;
894
895   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
896                stopwatch.RealTime(),stopwatch.CpuTime()));
897
898   return kTRUE;
899 }
900
901 //_____________________________________________________________________________
902 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
903 {
904 // run the HLT barrel tracking
905
906   TStopwatch stopwatch;
907   stopwatch.Start();
908
909   if (!fRunLoader) {
910     AliError("Missing runLoader!");
911     return kFALSE;
912   }
913
914   AliInfo("running HLT tracking");
915
916   // Get a pointer to the HLT reconstructor
917   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
918   if (!reconstructor) return kFALSE;
919
920   // TPC + ITS
921   for (Int_t iDet = 1; iDet >= 0; iDet--) {
922     TString detName = fgkDetectorName[iDet];
923     AliDebug(1, Form("%s HLT tracking", detName.Data()));
924     reconstructor->SetOption(detName.Data());
925     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
926     if (!tracker) {
927       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
928       if (fStopOnError) return kFALSE;
929       continue;
930     }
931     Double_t vtxPos[3];
932     Double_t vtxErr[3]={0.005,0.005,0.010};
933     const AliESDVertex *vertex = esd->GetVertex();
934     vertex->GetXYZ(vtxPos);
935     tracker->SetVertex(vtxPos,vtxErr);
936     if(iDet != 1) {
937       fLoader[iDet]->LoadRecPoints("read");
938       TTree* tree = fLoader[iDet]->TreeR();
939       if (!tree) {
940         AliError(Form("Can't get the %s cluster tree", detName.Data()));
941         return kFALSE;
942       }
943       tracker->LoadClusters(tree);
944     }
945     if (tracker->Clusters2Tracks(esd) != 0) {
946       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
947       return kFALSE;
948     }
949     if(iDet != 1) {
950       tracker->UnloadClusters();
951     }
952     delete tracker;
953   }
954
955   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
956                stopwatch.RealTime(),stopwatch.CpuTime()));
957
958   return kTRUE;
959 }
960
961 //_____________________________________________________________________________
962 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
963 {
964 // run the barrel tracking
965
966   TStopwatch stopwatch;
967   stopwatch.Start();
968
969   AliInfo("running tracking");
970
971   // pass 1: TPC + ITS inwards
972   for (Int_t iDet = 1; iDet >= 0; iDet--) {
973     if (!fTracker[iDet]) continue;
974     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
975
976     // load clusters
977     fLoader[iDet]->LoadRecPoints("read");
978     TTree* tree = fLoader[iDet]->TreeR();
979     if (!tree) {
980       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
981       return kFALSE;
982     }
983     fTracker[iDet]->LoadClusters(tree);
984
985     // run tracking
986     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
987       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
988       return kFALSE;
989     }
990     if (fCheckPointLevel > 1) {
991       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
992     }
993     // preliminary PID in TPC needed by the ITS tracker
994     if (iDet == 1) {
995       GetReconstructor(1)->FillESD(fRunLoader, esd);
996       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
997       AliESDpid::MakePID(esd);
998     }
999   }
1000
1001   // pass 2: ALL backwards
1002   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1003     if (!fTracker[iDet]) continue;
1004     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1005
1006     // load clusters
1007     if (iDet > 1) {     // all except ITS, TPC
1008       TTree* tree = NULL;
1009       fLoader[iDet]->LoadRecPoints("read");
1010       tree = fLoader[iDet]->TreeR();
1011       if (!tree) {
1012         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1013         return kFALSE;
1014       }
1015       fTracker[iDet]->LoadClusters(tree);
1016     }
1017
1018     // run tracking
1019     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1020       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1021       return kFALSE;
1022     }
1023     if (fCheckPointLevel > 1) {
1024       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1025     }
1026
1027     // unload clusters
1028     if (iDet > 2) {     // all except ITS, TPC, TRD
1029       fTracker[iDet]->UnloadClusters();
1030       fLoader[iDet]->UnloadRecPoints();
1031     }
1032     // updated PID in TPC needed by the ITS tracker -MI
1033     if (iDet == 1) {
1034       GetReconstructor(1)->FillESD(fRunLoader, esd);
1035       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1036       AliESDpid::MakePID(esd);
1037     }
1038   }
1039
1040   // write space-points to the ESD in case alignment data output
1041   // is switched on
1042   if (fWriteAlignmentData)
1043     WriteAlignmentData(esd);
1044
1045   // pass 3: TRD + TPC + ITS refit inwards
1046   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1047     if (!fTracker[iDet]) continue;
1048     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1049
1050     // run tracking
1051     if (fTracker[iDet]->RefitInward(esd) != 0) {
1052       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1053       return kFALSE;
1054     }
1055     if (fCheckPointLevel > 1) {
1056       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1057     }
1058
1059     // unload clusters
1060     fTracker[iDet]->UnloadClusters();
1061     fLoader[iDet]->UnloadRecPoints();
1062   }
1063   //
1064   // Propagate track to the vertex - if not done by ITS
1065   //
1066   Int_t ntracks = esd->GetNumberOfTracks();
1067   for (Int_t itrack=0; itrack<ntracks; itrack++){
1068     const Double_t kRadius  = 3;   // beam pipe radius
1069     const Double_t kMaxStep = 5;   // max step
1070     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1071     Double_t       fieldZ   = AliTracker::GetBz();  //
1072     AliESDtrack * track = esd->GetTrack(itrack);
1073     if (!track) continue;
1074     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1075     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1076     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1077   }
1078   
1079   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1080                stopwatch.RealTime(),stopwatch.CpuTime()));
1081
1082   return kTRUE;
1083 }
1084
1085 //_____________________________________________________________________________
1086 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1087 {
1088 // fill the event summary data
1089
1090   TStopwatch stopwatch;
1091   stopwatch.Start();
1092   AliInfo("filling ESD");
1093
1094   TString detStr = detectors;
1095   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1096     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1097     AliReconstructor* reconstructor = GetReconstructor(iDet);
1098     if (!reconstructor) continue;
1099
1100     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1101       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1102       TTree* clustersTree = NULL;
1103       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1104         fLoader[iDet]->LoadRecPoints("read");
1105         clustersTree = fLoader[iDet]->TreeR();
1106         if (!clustersTree) {
1107           AliError(Form("Can't get the %s clusters tree", 
1108                         fgkDetectorName[iDet]));
1109           if (fStopOnError) return kFALSE;
1110         }
1111       }
1112       if (fRawReader && !reconstructor->HasDigitConversion()) {
1113         reconstructor->FillESD(fRawReader, clustersTree, esd);
1114       } else {
1115         TTree* digitsTree = NULL;
1116         if (fLoader[iDet]) {
1117           fLoader[iDet]->LoadDigits("read");
1118           digitsTree = fLoader[iDet]->TreeD();
1119           if (!digitsTree) {
1120             AliError(Form("Can't get the %s digits tree", 
1121                           fgkDetectorName[iDet]));
1122             if (fStopOnError) return kFALSE;
1123           }
1124         }
1125         reconstructor->FillESD(digitsTree, clustersTree, esd);
1126         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1127       }
1128       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1129         fLoader[iDet]->UnloadRecPoints();
1130       }
1131
1132       if (fRawReader) {
1133         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1134       } else {
1135         reconstructor->FillESD(fRunLoader, esd);
1136       }
1137       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1138     }
1139   }
1140
1141   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1142     AliError(Form("the following detectors were not found: %s", 
1143                   detStr.Data()));
1144     if (fStopOnError) return kFALSE;
1145   }
1146
1147   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1148                stopwatch.RealTime(),stopwatch.CpuTime()));
1149
1150   return kTRUE;
1151 }
1152
1153 //_____________________________________________________________________________
1154 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1155 {
1156   // Reads the trigger decision which is
1157   // stored in Trigger.root file and fills
1158   // the corresponding esd entries
1159
1160   AliInfo("Filling trigger information into the ESD");
1161
1162   if (fRawReader) {
1163     AliCTPRawStream input(fRawReader);
1164     if (!input.Next()) {
1165       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1166       return kFALSE;
1167     }
1168     esd->SetTriggerMask(input.GetClassMask());
1169     esd->SetTriggerCluster(input.GetClusterMask());
1170   }
1171   else {
1172     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1173     if (runloader) {
1174       if (!runloader->LoadTrigger()) {
1175         AliCentralTrigger *aCTP = runloader->GetTrigger();
1176         esd->SetTriggerMask(aCTP->GetClassMask());
1177         esd->SetTriggerCluster(aCTP->GetClusterMask());
1178       }
1179       else {
1180         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1181         return kFALSE;
1182       }
1183     }
1184     else {
1185       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1186       return kFALSE;
1187     }
1188   }
1189
1190   return kTRUE;
1191 }
1192
1193 //_____________________________________________________________________________
1194 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1195 {
1196 // check whether detName is contained in detectors
1197 // if yes, it is removed from detectors
1198
1199   // check if all detectors are selected
1200   if ((detectors.CompareTo("ALL") == 0) ||
1201       detectors.BeginsWith("ALL ") ||
1202       detectors.EndsWith(" ALL") ||
1203       detectors.Contains(" ALL ")) {
1204     detectors = "ALL";
1205     return kTRUE;
1206   }
1207
1208   // search for the given detector
1209   Bool_t result = kFALSE;
1210   if ((detectors.CompareTo(detName) == 0) ||
1211       detectors.BeginsWith(detName+" ") ||
1212       detectors.EndsWith(" "+detName) ||
1213       detectors.Contains(" "+detName+" ")) {
1214     detectors.ReplaceAll(detName, "");
1215     result = kTRUE;
1216   }
1217
1218   // clean up the detectors string
1219   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1220   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1221   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1222
1223   return result;
1224 }
1225
1226 //_____________________________________________________________________________
1227 Bool_t AliReconstruction::InitRunLoader()
1228 {
1229 // get or create the run loader
1230
1231   if (gAlice) delete gAlice;
1232   gAlice = NULL;
1233
1234   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1235     // load all base libraries to get the loader classes
1236     TString libs = gSystem->GetLibraries();
1237     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1238       TString detName = fgkDetectorName[iDet];
1239       if (detName == "HLT") continue;
1240       if (libs.Contains("lib" + detName + "base.so")) continue;
1241       gSystem->Load("lib" + detName + "base.so");
1242     }
1243     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1244     if (!fRunLoader) {
1245       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1246       CleanUp();
1247       return kFALSE;
1248     }
1249     fRunLoader->CdGAFile();
1250     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1251       if (fRunLoader->LoadgAlice() == 0) {
1252         gAlice = fRunLoader->GetAliRun();
1253         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1254       }
1255     }
1256     if (!gAlice && !fRawReader) {
1257       AliError(Form("no gAlice object found in file %s",
1258                     fGAliceFileName.Data()));
1259       CleanUp();
1260       return kFALSE;
1261     }
1262
1263   } else {               // galice.root does not exist
1264     if (!fRawReader) {
1265       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1266       CleanUp();
1267       return kFALSE;
1268     }
1269     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1270                                     AliConfig::GetDefaultEventFolderName(),
1271                                     "recreate");
1272     if (!fRunLoader) {
1273       AliError(Form("could not create run loader in file %s", 
1274                     fGAliceFileName.Data()));
1275       CleanUp();
1276       return kFALSE;
1277     }
1278     fRunLoader->MakeTree("E");
1279     Int_t iEvent = 0;
1280     while (fRawReader->NextEvent()) {
1281       fRunLoader->SetEventNumber(iEvent);
1282       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
1283                                      iEvent, iEvent);
1284       fRunLoader->MakeTree("H");
1285       fRunLoader->TreeE()->Fill();
1286       iEvent++;
1287     }
1288     fRawReader->RewindEvents();
1289     fRunLoader->WriteHeader("OVERWRITE");
1290     fRunLoader->CdGAFile();
1291     fRunLoader->Write(0, TObject::kOverwrite);
1292 //    AliTracker::SetFieldMap(???);
1293   }
1294
1295   return kTRUE;
1296 }
1297
1298 //_____________________________________________________________________________
1299 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1300 {
1301 // get the reconstructor object and the loader for a detector
1302
1303   if (fReconstructor[iDet]) return fReconstructor[iDet];
1304
1305   // load the reconstructor object
1306   TPluginManager* pluginManager = gROOT->GetPluginManager();
1307   TString detName = fgkDetectorName[iDet];
1308   TString recName = "Ali" + detName + "Reconstructor";
1309   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1310
1311   if (detName == "HLT") {
1312     if (!gROOT->GetClass("AliLevel3")) {
1313       gSystem->Load("libAliL3Src.so");
1314       gSystem->Load("libAliL3Misc.so");
1315       gSystem->Load("libAliL3Hough.so");
1316       gSystem->Load("libAliL3Comp.so");
1317     }
1318   }
1319
1320   AliReconstructor* reconstructor = NULL;
1321   // first check if a plugin is defined for the reconstructor
1322   TPluginHandler* pluginHandler = 
1323     pluginManager->FindHandler("AliReconstructor", detName);
1324   // if not, add a plugin for it
1325   if (!pluginHandler) {
1326     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1327     TString libs = gSystem->GetLibraries();
1328     if (libs.Contains("lib" + detName + "base.so") ||
1329         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1330       pluginManager->AddHandler("AliReconstructor", detName, 
1331                                 recName, detName + "rec", recName + "()");
1332     } else {
1333       pluginManager->AddHandler("AliReconstructor", detName, 
1334                                 recName, detName, recName + "()");
1335     }
1336     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1337   }
1338   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1339     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1340   }
1341   if (reconstructor) {
1342     TObject* obj = fOptions.FindObject(detName.Data());
1343     if (obj) reconstructor->SetOption(obj->GetTitle());
1344     reconstructor->Init(fRunLoader);
1345     fReconstructor[iDet] = reconstructor;
1346   }
1347
1348   // get or create the loader
1349   if (detName != "HLT") {
1350     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1351     if (!fLoader[iDet]) {
1352       AliConfig::Instance()
1353         ->CreateDetectorFolders(fRunLoader->GetEventFolder(), 
1354                                 detName, detName);
1355       // first check if a plugin is defined for the loader
1356       TPluginHandler* pluginHandler = 
1357         pluginManager->FindHandler("AliLoader", detName);
1358       // if not, add a plugin for it
1359       if (!pluginHandler) {
1360         TString loaderName = "Ali" + detName + "Loader";
1361         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1362         pluginManager->AddHandler("AliLoader", detName, 
1363                                   loaderName, detName + "base", 
1364                                   loaderName + "(const char*, TFolder*)");
1365         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1366       }
1367       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1368         fLoader[iDet] = 
1369           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), 
1370                                                  fRunLoader->GetEventFolder());
1371       }
1372       if (!fLoader[iDet]) {   // use default loader
1373         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1374       }
1375       if (!fLoader[iDet]) {
1376         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1377         if (fStopOnError) return NULL;
1378       } else {
1379         fRunLoader->AddLoader(fLoader[iDet]);
1380         fRunLoader->CdGAFile();
1381         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1382         fRunLoader->Write(0, TObject::kOverwrite);
1383       }
1384     }
1385   }
1386       
1387   return reconstructor;
1388 }
1389
1390 //_____________________________________________________________________________
1391 Bool_t AliReconstruction::CreateVertexer()
1392 {
1393 // create the vertexer
1394
1395   fVertexer = NULL;
1396   AliReconstructor* itsReconstructor = GetReconstructor(0);
1397   if (itsReconstructor) {
1398     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1399   }
1400   if (!fVertexer) {
1401     AliWarning("couldn't create a vertexer for ITS");
1402     if (fStopOnError) return kFALSE;
1403   }
1404
1405   return kTRUE;
1406 }
1407
1408 //_____________________________________________________________________________
1409 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1410 {
1411 // create the trackers
1412
1413   TString detStr = detectors;
1414   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1415     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1416     AliReconstructor* reconstructor = GetReconstructor(iDet);
1417     if (!reconstructor) continue;
1418     TString detName = fgkDetectorName[iDet];
1419     if (detName == "HLT") {
1420       fRunHLTTracking = kTRUE;
1421       continue;
1422     }
1423
1424     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1425     if (!fTracker[iDet] && (iDet < 7)) {
1426       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1427       if (fStopOnError) return kFALSE;
1428     }
1429   }
1430
1431   return kTRUE;
1432 }
1433
1434 //_____________________________________________________________________________
1435 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1436 {
1437 // delete trackers and the run loader and close and delete the file
1438
1439   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1440     delete fReconstructor[iDet];
1441     fReconstructor[iDet] = NULL;
1442     fLoader[iDet] = NULL;
1443     delete fTracker[iDet];
1444     fTracker[iDet] = NULL;
1445   }
1446   delete fVertexer;
1447   fVertexer = NULL;
1448
1449   delete fRunLoader;
1450   fRunLoader = NULL;
1451   delete fRawReader;
1452   fRawReader = NULL;
1453
1454   if (file) {
1455     file->Close();
1456     delete file;
1457   }
1458
1459   if (fileOld) {
1460     fileOld->Close();
1461     delete fileOld;
1462     gSystem->Unlink("AliESDs.old.root");
1463   }
1464 }
1465
1466
1467 //_____________________________________________________________________________
1468 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1469 {
1470 // read the ESD event from a file
1471
1472   if (!esd) return kFALSE;
1473   char fileName[256];
1474   sprintf(fileName, "ESD_%d.%d_%s.root", 
1475           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1476   if (gSystem->AccessPathName(fileName)) return kFALSE;
1477
1478   AliInfo(Form("reading ESD from file %s", fileName));
1479   AliDebug(1, Form("reading ESD from file %s", fileName));
1480   TFile* file = TFile::Open(fileName);
1481   if (!file || !file->IsOpen()) {
1482     AliError(Form("opening %s failed", fileName));
1483     delete file;
1484     return kFALSE;
1485   }
1486
1487   gROOT->cd();
1488   delete esd;
1489   esd = (AliESD*) file->Get("ESD");
1490   file->Close();
1491   delete file;
1492   return kTRUE;
1493 }
1494
1495 //_____________________________________________________________________________
1496 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1497 {
1498 // write the ESD event to a file
1499
1500   if (!esd) return;
1501   char fileName[256];
1502   sprintf(fileName, "ESD_%d.%d_%s.root", 
1503           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1504
1505   AliDebug(1, Form("writing ESD to file %s", fileName));
1506   TFile* file = TFile::Open(fileName, "recreate");
1507   if (!file || !file->IsOpen()) {
1508     AliError(Form("opening %s failed", fileName));
1509   } else {
1510     esd->Write("ESD");
1511     file->Close();
1512   }
1513   delete file;
1514 }
1515
1516
1517
1518
1519 //_____________________________________________________________________________
1520 void AliReconstruction::CreateTag(TFile* file)
1521 {
1522   /////////////
1523   //muon code//
1524   ////////////
1525   Double_t fMUONMASS = 0.105658369;
1526   //Variables
1527   Double_t fX,fY,fZ ;
1528   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1529   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1530   Int_t fCharge;
1531   TLorentzVector fEPvector;
1532
1533   Float_t fZVertexCut = 10.0; 
1534   Float_t fRhoVertexCut = 2.0; 
1535
1536   Float_t fLowPtCut = 1.0;
1537   Float_t fHighPtCut = 3.0;
1538   Float_t fVeryHighPtCut = 10.0;
1539   ////////////
1540
1541   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1542
1543   // Creates the tags for all the events in a given ESD file
1544   Int_t ntrack;
1545   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1546   Int_t nPos, nNeg, nNeutr;
1547   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1548   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1549   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1550   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1551   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1552   Int_t fVertexflag;
1553   TString fVertexName("default");
1554
1555   AliRunTag *tag = new AliRunTag();
1556   AliEventTag *evTag = new AliEventTag();
1557   TTree ttag("T","A Tree with event tags");
1558   TBranch * btag = ttag.Branch("AliTAG", &tag);
1559   btag->SetCompressionLevel(9);
1560   
1561   AliInfo(Form("Creating the tags......."));    
1562   
1563   if (!file || !file->IsOpen()) {
1564     AliError(Form("opening failed"));
1565     delete file;
1566     return ;
1567   }  
1568   Int_t firstEvent = 0,lastEvent = 0;
1569   TTree *t = (TTree*) file->Get("esdTree");
1570   TBranch * b = t->GetBranch("ESD");
1571   AliESD *esd = 0;
1572   b->SetAddress(&esd);
1573   
1574   tag->SetRunId(esd->GetRunNumber());
1575   
1576   Int_t iNumberOfEvents = b->GetEntries();
1577   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1578     ntrack = 0;
1579     nPos = 0;
1580     nNeg = 0;
1581     nNeutr =0;
1582     nK0s = 0;
1583     nNeutrons = 0;
1584     nPi0s = 0;
1585     nGamas = 0;
1586     nProtons = 0;
1587     nKaons = 0;
1588     nPions = 0;
1589     nMuons = 0;
1590     nElectrons = 0;       
1591     nCh1GeV = 0;
1592     nCh3GeV = 0;
1593     nCh10GeV = 0;
1594     nMu1GeV = 0;
1595     nMu3GeV = 0;
1596     nMu10GeV = 0;
1597     nEl1GeV = 0;
1598     nEl3GeV = 0;
1599     nEl10GeV = 0;
1600     maxPt = .0;
1601     meanPt = .0;
1602     totalP = .0;
1603     fVertexflag = 0;
1604
1605     b->GetEntry(iEventNumber);
1606     const AliESDVertex * vertexIn = esd->GetVertex();
1607     if (!vertexIn) AliError("ESD has not defined vertex.");
1608     if (vertexIn) fVertexName = vertexIn->GetName();
1609     if(fVertexName != "default") fVertexflag = 1;
1610     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1611       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1612       UInt_t status = esdTrack->GetStatus();
1613       
1614       //select only tracks with ITS refit
1615       if ((status&AliESDtrack::kITSrefit)==0) continue;
1616       //select only tracks with TPC refit
1617       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1618       
1619       //select only tracks with the "combined PID"
1620       if ((status&AliESDtrack::kESDpid)==0) continue;
1621       Double_t p[3];
1622       esdTrack->GetPxPyPz(p);
1623       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1624       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1625       totalP += momentum;
1626       meanPt += fPt;
1627       if(fPt > maxPt) maxPt = fPt;
1628       
1629       if(esdTrack->GetSign() > 0) {
1630         nPos++;
1631         if(fPt > fLowPtCut) nCh1GeV++;
1632         if(fPt > fHighPtCut) nCh3GeV++;
1633         if(fPt > fVeryHighPtCut) nCh10GeV++;
1634       }
1635       if(esdTrack->GetSign() < 0) {
1636         nNeg++;
1637         if(fPt > fLowPtCut) nCh1GeV++;
1638         if(fPt > fHighPtCut) nCh3GeV++;
1639         if(fPt > fVeryHighPtCut) nCh10GeV++;
1640       }
1641       if(esdTrack->GetSign() == 0) nNeutr++;
1642       
1643       //PID
1644       Double_t prob[5];
1645       esdTrack->GetESDpid(prob);
1646       
1647       Double_t rcc = 0.0;
1648       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1649       if(rcc == 0.0) continue;
1650       //Bayes' formula
1651       Double_t w[5];
1652       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1653       
1654       //protons
1655       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1656       //kaons
1657       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1658       //pions
1659       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
1660       //electrons
1661       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1662         nElectrons++;
1663         if(fPt > fLowPtCut) nEl1GeV++;
1664         if(fPt > fHighPtCut) nEl3GeV++;
1665         if(fPt > fVeryHighPtCut) nEl10GeV++;
1666       }   
1667       ntrack++;
1668     }//track loop
1669     
1670     /////////////
1671     //muon code//
1672     ////////////
1673     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1674     // loop over all reconstructed tracks (also first track of combination)
1675     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1676       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1677       if (muonTrack == 0x0) continue;
1678       
1679       // Coordinates at vertex
1680       fZ = muonTrack->GetZ(); 
1681       fY = muonTrack->GetBendingCoor();
1682       fX = muonTrack->GetNonBendingCoor(); 
1683       
1684       fThetaX = muonTrack->GetThetaX();
1685       fThetaY = muonTrack->GetThetaY();
1686       
1687       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1688       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1689       fPxRec = fPzRec * TMath::Tan(fThetaX);
1690       fPyRec = fPzRec * TMath::Tan(fThetaY);
1691       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1692       
1693       //ChiSquare of the track if needed
1694       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1695       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1696       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1697       
1698       // total number of muons inside a vertex cut 
1699       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1700         nMuons++;
1701         if(fEPvector.Pt() > fLowPtCut) {
1702           nMu1GeV++; 
1703           if(fEPvector.Pt() > fHighPtCut) {
1704             nMu3GeV++; 
1705             if (fEPvector.Pt() > fVeryHighPtCut) {
1706               nMu10GeV++;
1707             }
1708           }
1709         }
1710       }
1711     }//muon track loop
1712     
1713     // Fill the event tags 
1714     if(ntrack != 0)
1715       meanPt = meanPt/ntrack;
1716     
1717     evTag->SetEventId(iEventNumber+1);
1718     if (vertexIn) {
1719       evTag->SetVertexX(vertexIn->GetXv());
1720       evTag->SetVertexY(vertexIn->GetYv());
1721       evTag->SetVertexZ(vertexIn->GetZv());
1722       evTag->SetVertexZError(vertexIn->GetZRes());
1723     }  
1724     evTag->SetVertexFlag(fVertexflag);
1725
1726     evTag->SetT0VertexZ(esd->GetT0zVertex());
1727     
1728     evTag->SetTrigger(esd->GetTriggerMask());
1729     
1730     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1731     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1732     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1733     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1734     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1735     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1736     
1737     
1738     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1739     evTag->SetNumOfPosTracks(nPos);
1740     evTag->SetNumOfNegTracks(nNeg);
1741     evTag->SetNumOfNeutrTracks(nNeutr);
1742     
1743     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1744     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1745     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1746     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1747     
1748     evTag->SetNumOfProtons(nProtons);
1749     evTag->SetNumOfKaons(nKaons);
1750     evTag->SetNumOfPions(nPions);
1751     evTag->SetNumOfMuons(nMuons);
1752     evTag->SetNumOfElectrons(nElectrons);
1753     evTag->SetNumOfPhotons(nGamas);
1754     evTag->SetNumOfPi0s(nPi0s);
1755     evTag->SetNumOfNeutrons(nNeutrons);
1756     evTag->SetNumOfKaon0s(nK0s);
1757     
1758     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1759     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1760     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1761     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1762     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1763     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1764     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1765     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1766     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1767     
1768     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1769     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1770     
1771     evTag->SetTotalMomentum(totalP);
1772     evTag->SetMeanPt(meanPt);
1773     evTag->SetMaxPt(maxPt);
1774     
1775     tag->AddEventTag(*evTag);
1776   }
1777   lastEvent = iNumberOfEvents;
1778         
1779   ttag.Fill();
1780   tag->Clear();
1781
1782   char fileName[256];
1783   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root", 
1784           tag->GetRunId(),firstEvent,lastEvent );
1785   AliInfo(Form("writing tags to file %s", fileName));
1786   AliDebug(1, Form("writing tags to file %s", fileName));
1787  
1788   TFile* ftag = TFile::Open(fileName, "recreate");
1789   ftag->cd();
1790   ttag.Write();
1791   ftag->Close();
1792   file->cd();
1793   delete tag;
1794   delete evTag;
1795 }
1796
1797 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1798 {
1799   // Write space-points which are then used in the alignment procedures
1800   // For the moment only ITS, TRD and TPC
1801
1802   // Load TOF clusters
1803   if (fTracker[3]){
1804     fLoader[3]->LoadRecPoints("read");
1805     TTree* tree = fLoader[3]->TreeR();
1806     if (!tree) {
1807       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1808       return;
1809     }
1810     fTracker[3]->LoadClusters(tree);
1811   }
1812   Int_t ntracks = esd->GetNumberOfTracks();
1813   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1814     {
1815       AliESDtrack *track = esd->GetTrack(itrack);
1816       Int_t nsp = 0;
1817       Int_t idx[200];
1818       for (Int_t iDet = 3; iDet >= 0; iDet--)
1819         nsp += track->GetNcls(iDet);
1820       if (nsp) {
1821         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1822         track->SetTrackPointArray(sp);
1823         Int_t isptrack = 0;
1824         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1825           AliTracker *tracker = fTracker[iDet];
1826           if (!tracker) continue;
1827           Int_t nspdet = track->GetNcls(iDet);
1828           if (nspdet <= 0) continue;
1829           track->GetClusters(iDet,idx);
1830           AliTrackPoint p;
1831           Int_t isp = 0;
1832           Int_t isp2 = 0;
1833           while (isp < nspdet) {
1834             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1835             const Int_t kNTPCmax = 159;
1836             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
1837             if (!isvalid) continue;
1838             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1839           }
1840         }       
1841       }
1842     }
1843   if (fTracker[3]){
1844     fTracker[3]->UnloadClusters();
1845     fLoader[3]->UnloadRecPoints();
1846   }
1847 }