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