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