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->GetTopNode()->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
 638     // Set magnetic field from the tracker
 639     esd->SetMagneticField(AliTracker::GetBz());
 640     hltesd->SetMagneticField(AliTracker::GetBz());
 641
 642     // vertex finder
 643     if (fRunVertexFinder) {
 644       if (!ReadESD(esd, "vertex")) {
 645         if (!RunVertexFinder(esd)) {
 646           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 647         }
 648         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 649       }
 650     }
 651
 652     // HLT tracking
 653     if (!fRunTracking.IsNull()) {
 654       if (fRunHLTTracking) {
 655         hltesd->SetVertex(esd->GetVertex());
 656         if (!RunHLTTracking(hltesd)) {
 657           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 658         }
 659       }
 660     }
 661
 662     // barrel tracking
 663     if (!fRunTracking.IsNull()) {
 664       if (!ReadESD(esd, "tracking")) {
 665         if (!RunTracking(esd)) {
 666           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 667         }
 668         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 669       }
 670     }
 671
 672     // fill ESD
 673     if (!fFillESD.IsNull()) {
 674       if (!FillESD(esd, fFillESD)) {
 675         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 676       }
 677     }
 678
 679     // combined PID
 680     AliESDpid::MakePID(esd);
 681     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 682
 683     if (fFillTriggerESD) {
 684       if (!ReadESD(esd, "trigger")) {
 685         if (!FillTriggerESD(esd)) {
 686           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 687         }
 688         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
 689       }
 690     }
 691
 692     esd->SetPrimaryVertex(tVertexer.FindPrimaryVertex(esd));
 693
 694     // write ESD
 695     tree->Fill();
 696     // write HLT ESD
 697     hlttree->Fill();
 698
 699     // write ESD friend
 700     if (fWriteESDfriend) {
 701        esdf=new AliESDfriend();
 702        esd->GetESDfriend(esdf);
 703        treef->Fill();
 704     }
 705
 706     if (fCheckPointLevel > 0)  WriteESD(esd, "final");
 707
 708     delete esd; delete hltesd;
 709     esd = NULL; hltesd = NULL;
 710   }
 711
 712   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
 713                stopwatch.RealTime(),stopwatch.CpuTime()));
 714
 715   file->cd();
 716   tree->Write();
 717   hlttree->Write();
 718
 719   if (fWriteESDfriend) {
 720      filef->cd();
 721      treef->Write(); delete treef; filef->Close(); delete filef;
 722   }
 723
 724   // Create tags for the events in the ESD tree (the ESD tree is always present)
 725   // In case of empty events the tags will contain dummy values
 726   CreateTag(file);
 727   CleanUp(file, fileOld);
 728
 729   return kTRUE;
 730 }
 731
 732
 733 //_____________________________________________________________________________
 734 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 735 {
 736 // run the local reconstruction
 737
 738   TStopwatch stopwatch;
 739   stopwatch.Start();
 740
 741   TString detStr = detectors;
 742   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 743     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 744     AliReconstructor* reconstructor = GetReconstructor(iDet);
 745     if (!reconstructor) continue;
 746     if (reconstructor->HasLocalReconstruction()) continue;
 747
 748     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 749     TStopwatch stopwatchDet;
 750     stopwatchDet.Start();
 751     if (fRawReader) {
 752       fRawReader->RewindEvents();
 753       reconstructor->Reconstruct(fRunLoader, fRawReader);
 754     } else {
 755       reconstructor->Reconstruct(fRunLoader);
 756     }
 757     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 758                  fgkDetectorName[iDet],
 759                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 760   }
 761
 762   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 763     AliError(Form("the following detectors were not found: %s",
 764                   detStr.Data()));
 765     if (fStopOnError) return kFALSE;
 766   }
 767
 768   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 769                stopwatch.RealTime(),stopwatch.CpuTime()));
 770
 771   return kTRUE;
 772 }
 773
 774 //_____________________________________________________________________________
 775 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 776 {
 777 // run the local reconstruction
 778
 779   TStopwatch stopwatch;
 780   stopwatch.Start();
 781
 782   TString detStr = detectors;
 783   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 784     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 785     AliReconstructor* reconstructor = GetReconstructor(iDet);
 786     if (!reconstructor) continue;
 787     AliLoader* loader = fLoader[iDet];
 788
 789     // conversion of digits
 790     if (fRawReader && reconstructor->HasDigitConversion()) {
 791       AliInfo(Form("converting raw data digits into root objects for %s",
 792                    fgkDetectorName[iDet]));
 793       TStopwatch stopwatchDet;
 794       stopwatchDet.Start();
 795       loader->LoadDigits("update");
 796       loader->CleanDigits();
 797       loader->MakeDigitsContainer();
 798       TTree* digitsTree = loader->TreeD();
 799       reconstructor->ConvertDigits(fRawReader, digitsTree);
 800       loader->WriteDigits("OVERWRITE");
 801       loader->UnloadDigits();
 802       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 803                    fgkDetectorName[iDet],
 804                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 805     }
 806
 807     // local reconstruction
 808     if (!reconstructor->HasLocalReconstruction()) continue;
 809     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 810     TStopwatch stopwatchDet;
 811     stopwatchDet.Start();
 812     loader->LoadRecPoints("update");
 813     loader->CleanRecPoints();
 814     loader->MakeRecPointsContainer();
 815     TTree* clustersTree = loader->TreeR();
 816     if (fRawReader && !reconstructor->HasDigitConversion()) {
 817       reconstructor->Reconstruct(fRawReader, clustersTree);
 818     } else {
 819       loader->LoadDigits("read");
 820       TTree* digitsTree = loader->TreeD();
 821       if (!digitsTree) {
 822         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 823         if (fStopOnError) return kFALSE;
 824       } else {
 825         reconstructor->Reconstruct(digitsTree, clustersTree);
 826       }
 827       loader->UnloadDigits();
 828     }
 829     loader->WriteRecPoints("OVERWRITE");
 830     loader->UnloadRecPoints();
 831     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
 832                     fgkDetectorName[iDet],
 833                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 834   }
 835
 836   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 837     AliError(Form("the following detectors were not found: %s",
 838                   detStr.Data()));
 839     if (fStopOnError) return kFALSE;
 840   }
 841
 842   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 843                stopwatch.RealTime(),stopwatch.CpuTime()));
 844
 845   return kTRUE;
 846 }
 847
 848 //_____________________________________________________________________________
 849 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 850 {
 851 // run the barrel tracking
 852
 853   TStopwatch stopwatch;
 854   stopwatch.Start();
 855
 856   AliESDVertex* vertex = NULL;
 857   Double_t vtxPos[3] = {0, 0, 0};
 858   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 859   TArrayF mcVertex(3);
 860   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 861     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 862     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 863   }
 864
 865   if (fVertexer) {
 866     AliInfo("running the ITS vertex finder");
 867     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 868     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 869     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 870     if(!vertex){
 871       AliWarning("Vertex not found");
 872       vertex = new AliESDVertex();
 873       vertex->SetName("default");
 874     }
 875     else {
 876       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 877       vertex->SetName("reconstructed");
 878     }
 879
 880   } else {
 881     AliInfo("getting the primary vertex from MC");
 882     vertex = new AliESDVertex(vtxPos, vtxErr);
 883   }
 884
 885   if (vertex) {
 886     vertex->GetXYZ(vtxPos);
 887     vertex->GetSigmaXYZ(vtxErr);
 888   } else {
 889     AliWarning("no vertex reconstructed");
 890     vertex = new AliESDVertex(vtxPos, vtxErr);
 891   }
 892   esd->SetVertex(vertex);
 893   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 894     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
 895   }
 896   delete vertex;
 897
 898   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 899                stopwatch.RealTime(),stopwatch.CpuTime()));
 900
 901   return kTRUE;
 902 }
 903
 904 //_____________________________________________________________________________
 905 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
 906 {
 907 // run the HLT barrel tracking
 908
 909   TStopwatch stopwatch;
 910   stopwatch.Start();
 911
 912   if (!fRunLoader) {
 913     AliError("Missing runLoader!");
 914     return kFALSE;
 915   }
 916
 917   AliInfo("running HLT tracking");
 918
 919   // Get a pointer to the HLT reconstructor
 920   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
 921   if (!reconstructor) return kFALSE;
 922
 923   // TPC + ITS
 924   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 925     TString detName = fgkDetectorName[iDet];
 926     AliDebug(1, Form("%s HLT tracking", detName.Data()));
 927     reconstructor->SetOption(detName.Data());
 928     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
 929     if (!tracker) {
 930       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
 931       if (fStopOnError) return kFALSE;
 932       continue;
 933     }
 934     Double_t vtxPos[3];
 935     Double_t vtxErr[3]={0.005,0.005,0.010};
 936     const AliESDVertex *vertex = esd->GetVertex();
 937     vertex->GetXYZ(vtxPos);
 938     tracker->SetVertex(vtxPos,vtxErr);
 939     if(iDet != 1) {
 940       fLoader[iDet]->LoadRecPoints("read");
 941       TTree* tree = fLoader[iDet]->TreeR();
 942       if (!tree) {
 943         AliError(Form("Can't get the %s cluster tree", detName.Data()));
 944         return kFALSE;
 945       }
 946       tracker->LoadClusters(tree);
 947     }
 948     if (tracker->Clusters2Tracks(esd) != 0) {
 949       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
 950       return kFALSE;
 951     }
 952     if(iDet != 1) {
 953       tracker->UnloadClusters();
 954     }
 955     delete tracker;
 956   }
 957
 958   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 959                stopwatch.RealTime(),stopwatch.CpuTime()));
 960
 961   return kTRUE;
 962 }
 963
 964 //_____________________________________________________________________________
 965 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
 966 {
 967 // run the barrel tracking
 968
 969   TStopwatch stopwatch;
 970   stopwatch.Start();
 971
 972   AliInfo("running tracking");
 973
 974   // pass 1: TPC + ITS inwards
 975   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 976     if (!fTracker[iDet]) continue;
 977     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
 978
 979     // load clusters
 980     fLoader[iDet]->LoadRecPoints("read");
 981     TTree* tree = fLoader[iDet]->TreeR();
 982     if (!tree) {
 983       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 984       return kFALSE;
 985     }
 986     fTracker[iDet]->LoadClusters(tree);
 987
 988     // run tracking
 989     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
 990       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
 991       return kFALSE;
 992     }
 993     if (fCheckPointLevel > 1) {
 994       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
 995     }
 996     // preliminary PID in TPC needed by the ITS tracker
 997     if (iDet == 1) {
 998       GetReconstructor(1)->FillESD(fRunLoader, esd);
 999       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1000       AliESDpid::MakePID(esd);
1001     }
1002   }
1003
1004   // pass 2: ALL backwards
1005   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1006     if (!fTracker[iDet]) continue;
1007     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1008
1009     // load clusters
1010     if (iDet > 1) {     // all except ITS, TPC
1011       TTree* tree = NULL;
1012       fLoader[iDet]->LoadRecPoints("read");
1013       tree = fLoader[iDet]->TreeR();
1014       if (!tree) {
1015         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1016         return kFALSE;
1017       }
1018       fTracker[iDet]->LoadClusters(tree);
1019     }
1020
1021     // run tracking
1022     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1023       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1024       return kFALSE;
1025     }
1026     if (fCheckPointLevel > 1) {
1027       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1028     }
1029
1030     // unload clusters
1031     if (iDet > 2) {     // all except ITS, TPC, TRD
1032       fTracker[iDet]->UnloadClusters();
1033       fLoader[iDet]->UnloadRecPoints();
1034     }
1035     // updated PID in TPC needed by the ITS tracker -MI
1036     if (iDet == 1) {
1037       GetReconstructor(1)->FillESD(fRunLoader, esd);
1038       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1039       AliESDpid::MakePID(esd);
1040     }
1041   }
1042
1043   // write space-points to the ESD in case alignment data output
1044   // is switched on
1045   if (fWriteAlignmentData)
1046     WriteAlignmentData(esd);
1047
1048   // pass 3: TRD + TPC + ITS refit inwards
1049   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1050     if (!fTracker[iDet]) continue;
1051     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1052
1053     // run tracking
1054     if (fTracker[iDet]->RefitInward(esd) != 0) {
1055       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1056       return kFALSE;
1057     }
1058     if (fCheckPointLevel > 1) {
1059       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1060     }
1061
1062     // unload clusters
1063     fTracker[iDet]->UnloadClusters();
1064     fLoader[iDet]->UnloadRecPoints();
1065   }
1066   //
1067   // Propagate track to the vertex - if not done by ITS
1068   //
1069   Int_t ntracks = esd->GetNumberOfTracks();
1070   for (Int_t itrack=0; itrack<ntracks; itrack++){
1071     const Double_t kRadius  = 3;   // beam pipe radius
1072     const Double_t kMaxStep = 5;   // max step
1073     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1074     Double_t       fieldZ   = AliTracker::GetBz();  //
1075     AliESDtrack * track = esd->GetTrack(itrack);
1076     if (!track) continue;
1077     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1078     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1079     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1080   }
1081
1082   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1083                stopwatch.RealTime(),stopwatch.CpuTime()));
1084
1085   return kTRUE;
1086 }
1087
1088 //_____________________________________________________________________________
1089 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1090 {
1091 // fill the event summary data
1092
1093   TStopwatch stopwatch;
1094   stopwatch.Start();
1095   AliInfo("filling ESD");
1096
1097   TString detStr = detectors;
1098   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1099     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1100     AliReconstructor* reconstructor = GetReconstructor(iDet);
1101     if (!reconstructor) continue;
1102
1103     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1104       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1105       TTree* clustersTree = NULL;
1106       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1107         fLoader[iDet]->LoadRecPoints("read");
1108         clustersTree = fLoader[iDet]->TreeR();
1109         if (!clustersTree) {
1110           AliError(Form("Can't get the %s clusters tree",
1111                         fgkDetectorName[iDet]));
1112           if (fStopOnError) return kFALSE;
1113         }
1114       }
1115       if (fRawReader && !reconstructor->HasDigitConversion()) {
1116         reconstructor->FillESD(fRawReader, clustersTree, esd);
1117       } else {
1118         TTree* digitsTree = NULL;
1119         if (fLoader[iDet]) {
1120           fLoader[iDet]->LoadDigits("read");
1121           digitsTree = fLoader[iDet]->TreeD();
1122           if (!digitsTree) {
1123             AliError(Form("Can't get the %s digits tree",
1124                           fgkDetectorName[iDet]));
1125             if (fStopOnError) return kFALSE;
1126           }
1127         }
1128         reconstructor->FillESD(digitsTree, clustersTree, esd);
1129         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1130       }
1131       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1132         fLoader[iDet]->UnloadRecPoints();
1133       }
1134
1135       if (fRawReader) {
1136         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1137       } else {
1138         reconstructor->FillESD(fRunLoader, esd);
1139       }
1140       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1141     }
1142   }
1143
1144   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1145     AliError(Form("the following detectors were not found: %s",
1146                   detStr.Data()));
1147     if (fStopOnError) return kFALSE;
1148   }
1149
1150   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1151                stopwatch.RealTime(),stopwatch.CpuTime()));
1152
1153   return kTRUE;
1154 }
1155
1156 //_____________________________________________________________________________
1157 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1158 {
1159   // Reads the trigger decision which is
1160   // stored in Trigger.root file and fills
1161   // the corresponding esd entries
1162
1163   AliInfo("Filling trigger information into the ESD");
1164
1165   if (fRawReader) {
1166     AliCTPRawStream input(fRawReader);
1167     if (!input.Next()) {
1168       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1169       return kFALSE;
1170     }
1171     esd->SetTriggerMask(input.GetClassMask());
1172     esd->SetTriggerCluster(input.GetClusterMask());
1173   }
1174   else {
1175     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1176     if (runloader) {
1177       if (!runloader->LoadTrigger()) {
1178         AliCentralTrigger *aCTP = runloader->GetTrigger();
1179         esd->SetTriggerMask(aCTP->GetClassMask());
1180         esd->SetTriggerCluster(aCTP->GetClusterMask());
1181       }
1182       else {
1183         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1184         return kFALSE;
1185       }
1186     }
1187     else {
1188       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1189       return kFALSE;
1190     }
1191   }
1192
1193   return kTRUE;
1194 }
1195
1196 //_____________________________________________________________________________
1197 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1198 {
1199 // check whether detName is contained in detectors
1200 // if yes, it is removed from detectors
1201
1202   // check if all detectors are selected
1203   if ((detectors.CompareTo("ALL") == 0) ||
1204       detectors.BeginsWith("ALL ") ||
1205       detectors.EndsWith(" ALL") ||
1206       detectors.Contains(" ALL ")) {
1207     detectors = "ALL";
1208     return kTRUE;
1209   }
1210
1211   // search for the given detector
1212   Bool_t result = kFALSE;
1213   if ((detectors.CompareTo(detName) == 0) ||
1214       detectors.BeginsWith(detName+" ") ||
1215       detectors.EndsWith(" "+detName) ||
1216       detectors.Contains(" "+detName+" ")) {
1217     detectors.ReplaceAll(detName, "");
1218     result = kTRUE;
1219   }
1220
1221   // clean up the detectors string
1222   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1223   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1224   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1225
1226   return result;
1227 }
1228
1229 //_____________________________________________________________________________
1230 Bool_t AliReconstruction::InitRunLoader()
1231 {
1232 // get or create the run loader
1233
1234   if (gAlice) delete gAlice;
1235   gAlice = NULL;
1236
1237   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1238     // load all base libraries to get the loader classes
1239     TString libs = gSystem->GetLibraries();
1240     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1241       TString detName = fgkDetectorName[iDet];
1242       if (detName == "HLT") continue;
1243       if (libs.Contains("lib" + detName + "base.so")) continue;
1244       gSystem->Load("lib" + detName + "base.so");
1245     }
1246     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1247     if (!fRunLoader) {
1248       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1249       CleanUp();
1250       return kFALSE;
1251     }
1252     fRunLoader->CdGAFile();
1253     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1254       if (fRunLoader->LoadgAlice() == 0) {
1255         gAlice = fRunLoader->GetAliRun();
1256         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1257       }
1258     }
1259     if (!gAlice && !fRawReader) {
1260       AliError(Form("no gAlice object found in file %s",
1261                     fGAliceFileName.Data()));
1262       CleanUp();
1263       return kFALSE;
1264     }
1265
1266   } else {               // galice.root does not exist
1267     if (!fRawReader) {
1268       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1269       CleanUp();
1270       return kFALSE;
1271     }
1272     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1273                                     AliConfig::GetDefaultEventFolderName(),
1274                                     "recreate");
1275     if (!fRunLoader) {
1276       AliError(Form("could not create run loader in file %s",
1277                     fGAliceFileName.Data()));
1278       CleanUp();
1279       return kFALSE;
1280     }
1281     fRunLoader->MakeTree("E");
1282     Int_t iEvent = 0;
1283     while (fRawReader->NextEvent()) {
1284       fRunLoader->SetEventNumber(iEvent);
1285       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1286                                      iEvent, iEvent);
1287       fRunLoader->MakeTree("H");
1288       fRunLoader->TreeE()->Fill();
1289       iEvent++;
1290     }
1291     fRawReader->RewindEvents();
1292     fRunLoader->WriteHeader("OVERWRITE");
1293     fRunLoader->CdGAFile();
1294     fRunLoader->Write(0, TObject::kOverwrite);
1295 //    AliTracker::SetFieldMap(???);
1296   }
1297
1298   return kTRUE;
1299 }
1300
1301 //_____________________________________________________________________________
1302 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1303 {
1304 // get the reconstructor object and the loader for a detector
1305
1306   if (fReconstructor[iDet]) return fReconstructor[iDet];
1307
1308   // load the reconstructor object
1309   TPluginManager* pluginManager = gROOT->GetPluginManager();
1310   TString detName = fgkDetectorName[iDet];
1311   TString recName = "Ali" + detName + "Reconstructor";
1312   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1313
1314   if (detName == "HLT") {
1315     if (!gROOT->GetClass("AliLevel3")) {
1316       gSystem->Load("libAliL3Src.so");
1317       gSystem->Load("libAliL3Misc.so");
1318       gSystem->Load("libAliL3Hough.so");
1319       gSystem->Load("libAliL3Comp.so");
1320     }
1321   }
1322
1323   AliReconstructor* reconstructor = NULL;
1324   // first check if a plugin is defined for the reconstructor
1325   TPluginHandler* pluginHandler =
1326     pluginManager->FindHandler("AliReconstructor", detName);
1327   // if not, add a plugin for it
1328   if (!pluginHandler) {
1329     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1330     TString libs = gSystem->GetLibraries();
1331     if (libs.Contains("lib" + detName + "base.so") ||
1332         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1333       pluginManager->AddHandler("AliReconstructor", detName,
1334                                 recName, detName + "rec", recName + "()");
1335     } else {
1336       pluginManager->AddHandler("AliReconstructor", detName,
1337                                 recName, detName, recName + "()");
1338     }
1339     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1340   }
1341   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1342     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1343   }
1344   if (reconstructor) {
1345     TObject* obj = fOptions.FindObject(detName.Data());
1346     if (obj) reconstructor->SetOption(obj->GetTitle());
1347     reconstructor->Init(fRunLoader);
1348     fReconstructor[iDet] = reconstructor;
1349   }
1350
1351   // get or create the loader
1352   if (detName != "HLT") {
1353     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1354     if (!fLoader[iDet]) {
1355       AliConfig::Instance()
1356         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1357                                 detName, detName);
1358       // first check if a plugin is defined for the loader
1359       TPluginHandler* pluginHandler =
1360         pluginManager->FindHandler("AliLoader", detName);
1361       // if not, add a plugin for it
1362       if (!pluginHandler) {
1363         TString loaderName = "Ali" + detName + "Loader";
1364         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1365         pluginManager->AddHandler("AliLoader", detName,
1366                                   loaderName, detName + "base",
1367                                   loaderName + "(const char*, TFolder*)");
1368         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1369       }
1370       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1371         fLoader[iDet] =
1372           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1373                                                  fRunLoader->GetEventFolder());
1374       }
1375       if (!fLoader[iDet]) {   // use default loader
1376         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1377       }
1378       if (!fLoader[iDet]) {
1379         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1380         if (fStopOnError) return NULL;
1381       } else {
1382         fRunLoader->AddLoader(fLoader[iDet]);
1383         fRunLoader->CdGAFile();
1384         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1385         fRunLoader->Write(0, TObject::kOverwrite);
1386       }
1387     }
1388   }
1389
1390   return reconstructor;
1391 }
1392
1393 //_____________________________________________________________________________
1394 Bool_t AliReconstruction::CreateVertexer()
1395 {
1396 // create the vertexer
1397
1398   fVertexer = NULL;
1399   AliReconstructor* itsReconstructor = GetReconstructor(0);
1400   if (itsReconstructor) {
1401     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1402   }
1403   if (!fVertexer) {
1404     AliWarning("couldn't create a vertexer for ITS");
1405     if (fStopOnError) return kFALSE;
1406   }
1407
1408   return kTRUE;
1409 }
1410
1411 //_____________________________________________________________________________
1412 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1413 {
1414 // create the trackers
1415
1416   TString detStr = detectors;
1417   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1418     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1419     AliReconstructor* reconstructor = GetReconstructor(iDet);
1420     if (!reconstructor) continue;
1421     TString detName = fgkDetectorName[iDet];
1422     if (detName == "HLT") {
1423       fRunHLTTracking = kTRUE;
1424       continue;
1425     }
1426
1427     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1428     if (!fTracker[iDet] && (iDet < 7)) {
1429       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1430       if (fStopOnError) return kFALSE;
1431     }
1432   }
1433
1434   return kTRUE;
1435 }
1436
1437 //_____________________________________________________________________________
1438 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1439 {
1440 // delete trackers and the run loader and close and delete the file
1441
1442   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1443     delete fReconstructor[iDet];
1444     fReconstructor[iDet] = NULL;
1445     fLoader[iDet] = NULL;
1446     delete fTracker[iDet];
1447     fTracker[iDet] = NULL;
1448   }
1449   delete fVertexer;
1450   fVertexer = NULL;
1451
1452   delete fRunLoader;
1453   fRunLoader = NULL;
1454   delete fRawReader;
1455   fRawReader = NULL;
1456
1457   if (file) {
1458     file->Close();
1459     delete file;
1460   }
1461
1462   if (fileOld) {
1463     fileOld->Close();
1464     delete fileOld;
1465     gSystem->Unlink("AliESDs.old.root");
1466   }
1467 }
1468
1469
1470 //_____________________________________________________________________________
1471 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1472 {
1473 // read the ESD event from a file
1474
1475   if (!esd) return kFALSE;
1476   char fileName[256];
1477   sprintf(fileName, "ESD_%d.%d_%s.root",
1478           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1479   if (gSystem->AccessPathName(fileName)) return kFALSE;
1480
1481   AliInfo(Form("reading ESD from file %s", fileName));
1482   AliDebug(1, Form("reading ESD from file %s", fileName));
1483   TFile* file = TFile::Open(fileName);
1484   if (!file || !file->IsOpen()) {
1485     AliError(Form("opening %s failed", fileName));
1486     delete file;
1487     return kFALSE;
1488   }
1489
1490   gROOT->cd();
1491   delete esd;
1492   esd = (AliESD*) file->Get("ESD");
1493   file->Close();
1494   delete file;
1495   return kTRUE;
1496 }
1497
1498 //_____________________________________________________________________________
1499 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1500 {
1501 // write the ESD event to a file
1502
1503   if (!esd) return;
1504   char fileName[256];
1505   sprintf(fileName, "ESD_%d.%d_%s.root",
1506           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1507
1508   AliDebug(1, Form("writing ESD to file %s", fileName));
1509   TFile* file = TFile::Open(fileName, "recreate");
1510   if (!file || !file->IsOpen()) {
1511     AliError(Form("opening %s failed", fileName));
1512   } else {
1513     esd->Write("ESD");
1514     file->Close();
1515   }
1516   delete file;
1517 }
1518
1519
1520
1521
1522 //_____________________________________________________________________________
1523 void AliReconstruction::CreateTag(TFile* file)
1524 {
1525   /////////////
1526   //muon code//
1527   ////////////
1528   Double_t fMUONMASS = 0.105658369;
1529   //Variables
1530   Double_t fX,fY,fZ ;
1531   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1532   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1533   Int_t fCharge;
1534   TLorentzVector fEPvector;
1535
1536   Float_t fZVertexCut = 10.0;
1537   Float_t fRhoVertexCut = 2.0;
1538
1539   Float_t fLowPtCut = 1.0;
1540   Float_t fHighPtCut = 3.0;
1541   Float_t fVeryHighPtCut = 10.0;
1542   ////////////
1543
1544   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1545
1546   // Creates the tags for all the events in a given ESD file
1547   Int_t ntrack;
1548   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1549   Int_t nPos, nNeg, nNeutr;
1550   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1551   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1552   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1553   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1554   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1555   Int_t fVertexflag;
1556   Int_t iRunNumber = 0;
1557   TString fVertexName("default");
1558
1559   AliRunTag *tag = new AliRunTag();
1560   AliEventTag *evTag = new AliEventTag();
1561   TTree ttag("T","A Tree with event tags");
1562   TBranch * btag = ttag.Branch("AliTAG", &tag);
1563   btag->SetCompressionLevel(9);
1564
1565   AliInfo(Form("Creating the tags......."));
1566
1567   if (!file || !file->IsOpen()) {
1568     AliError(Form("opening failed"));
1569     delete file;
1570     return ;
1571   }
1572   Int_t firstEvent = 0,lastEvent = 0;
1573   TTree *t = (TTree*) file->Get("esdTree");
1574   TBranch * b = t->GetBranch("ESD");
1575   AliESD *esd = 0;
1576   b->SetAddress(&esd);
1577
1578   b->GetEntry(0);
1579   Int_t iInitRunNumber = 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     iRunNumber = esd->GetRunNumber();
1612     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1613     const AliESDVertex * vertexIn = esd->GetVertex();
1614     if (!vertexIn) AliError("ESD has not defined vertex.");
1615     if (vertexIn) fVertexName = vertexIn->GetName();
1616     if(fVertexName != "default") fVertexflag = 1;
1617     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1618       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1619       UInt_t status = esdTrack->GetStatus();
1620
1621       //select only tracks with ITS refit
1622       if ((status&AliESDtrack::kITSrefit)==0) continue;
1623       //select only tracks with TPC refit
1624       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1625
1626       //select only tracks with the "combined PID"
1627       if ((status&AliESDtrack::kESDpid)==0) continue;
1628       Double_t p[3];
1629       esdTrack->GetPxPyPz(p);
1630       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1631       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1632       totalP += momentum;
1633       meanPt += fPt;
1634       if(fPt > maxPt) maxPt = fPt;
1635
1636       if(esdTrack->GetSign() > 0) {
1637         nPos++;
1638         if(fPt > fLowPtCut) nCh1GeV++;
1639         if(fPt > fHighPtCut) nCh3GeV++;
1640         if(fPt > fVeryHighPtCut) nCh10GeV++;
1641       }
1642       if(esdTrack->GetSign() < 0) {
1643         nNeg++;
1644         if(fPt > fLowPtCut) nCh1GeV++;
1645         if(fPt > fHighPtCut) nCh3GeV++;
1646         if(fPt > fVeryHighPtCut) nCh10GeV++;
1647       }
1648       if(esdTrack->GetSign() == 0) nNeutr++;
1649
1650       //PID
1651       Double_t prob[5];
1652       esdTrack->GetESDpid(prob);
1653
1654       Double_t rcc = 0.0;
1655       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1656       if(rcc == 0.0) continue;
1657       //Bayes' formula
1658       Double_t w[5];
1659       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1660
1661       //protons
1662       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1663       //kaons
1664       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1665       //pions
1666       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1667       //electrons
1668       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1669         nElectrons++;
1670         if(fPt > fLowPtCut) nEl1GeV++;
1671         if(fPt > fHighPtCut) nEl3GeV++;
1672         if(fPt > fVeryHighPtCut) nEl10GeV++;
1673       }
1674       ntrack++;
1675     }//track loop
1676
1677     /////////////
1678     //muon code//
1679     ////////////
1680     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1681     // loop over all reconstructed tracks (also first track of combination)
1682     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1683       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1684       if (muonTrack == 0x0) continue;
1685
1686       // Coordinates at vertex
1687       fZ = muonTrack->GetZ();
1688       fY = muonTrack->GetBendingCoor();
1689       fX = muonTrack->GetNonBendingCoor();
1690
1691       fThetaX = muonTrack->GetThetaX();
1692       fThetaY = muonTrack->GetThetaY();
1693
1694       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1695       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1696       fPxRec = fPzRec * TMath::Tan(fThetaX);
1697       fPyRec = fPzRec * TMath::Tan(fThetaY);
1698       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1699
1700       //ChiSquare of the track if needed
1701       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1702       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1703       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1704
1705       // total number of muons inside a vertex cut
1706       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1707         nMuons++;
1708         if(fEPvector.Pt() > fLowPtCut) {
1709           nMu1GeV++;
1710           if(fEPvector.Pt() > fHighPtCut) {
1711             nMu3GeV++;
1712             if (fEPvector.Pt() > fVeryHighPtCut) {
1713               nMu10GeV++;
1714             }
1715           }
1716         }
1717       }
1718     }//muon track loop
1719
1720     // Fill the event tags
1721     if(ntrack != 0)
1722       meanPt = meanPt/ntrack;
1723
1724     evTag->SetEventId(iEventNumber+1);
1725     if (vertexIn) {
1726       evTag->SetVertexX(vertexIn->GetXv());
1727       evTag->SetVertexY(vertexIn->GetYv());
1728       evTag->SetVertexZ(vertexIn->GetZv());
1729       evTag->SetVertexZError(vertexIn->GetZRes());
1730     }
1731     evTag->SetVertexFlag(fVertexflag);
1732
1733     evTag->SetT0VertexZ(esd->GetT0zVertex());
1734
1735     evTag->SetTriggerMask(esd->GetTriggerMask());
1736     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1737
1738     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1739     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1740     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1741     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1742     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1743     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1744
1745
1746     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1747     evTag->SetNumOfPosTracks(nPos);
1748     evTag->SetNumOfNegTracks(nNeg);
1749     evTag->SetNumOfNeutrTracks(nNeutr);
1750
1751     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1752     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1753     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1754     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1755
1756     evTag->SetNumOfProtons(nProtons);
1757     evTag->SetNumOfKaons(nKaons);
1758     evTag->SetNumOfPions(nPions);
1759     evTag->SetNumOfMuons(nMuons);
1760     evTag->SetNumOfElectrons(nElectrons);
1761     evTag->SetNumOfPhotons(nGamas);
1762     evTag->SetNumOfPi0s(nPi0s);
1763     evTag->SetNumOfNeutrons(nNeutrons);
1764     evTag->SetNumOfKaon0s(nK0s);
1765
1766     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1767     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1768     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1769     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1770     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1771     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1772     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1773     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1774     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1775
1776     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1777     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1778
1779     evTag->SetTotalMomentum(totalP);
1780     evTag->SetMeanPt(meanPt);
1781     evTag->SetMaxPt(maxPt);
1782
1783     tag->SetRunId(iInitRunNumber);
1784     tag->AddEventTag(*evTag);
1785   }
1786   lastEvent = iNumberOfEvents;
1787
1788   ttag.Fill();
1789   tag->Clear();
1790
1791   char fileName[256];
1792   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1793           tag->GetRunId(),firstEvent,lastEvent );
1794   AliInfo(Form("writing tags to file %s", fileName));
1795   AliDebug(1, Form("writing tags to file %s", fileName));
1796
1797   TFile* ftag = TFile::Open(fileName, "recreate");
1798   ftag->cd();
1799   ttag.Write();
1800   ftag->Close();
1801   file->cd();
1802   delete tag;
1803   delete evTag;
1804 }
1805
1806 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1807 {
1808   // Write space-points which are then used in the alignment procedures
1809   // For the moment only ITS, TRD and TPC
1810
1811   // Load TOF clusters
1812   if (fTracker[3]){
1813     fLoader[3]->LoadRecPoints("read");
1814     TTree* tree = fLoader[3]->TreeR();
1815     if (!tree) {
1816       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1817       return;
1818     }
1819     fTracker[3]->LoadClusters(tree);
1820   }
1821   Int_t ntracks = esd->GetNumberOfTracks();
1822   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1823     {
1824       AliESDtrack *track = esd->GetTrack(itrack);
1825       Int_t nsp = 0;
1826       Int_t idx[200];
1827       for (Int_t iDet = 3; iDet >= 0; iDet--)
1828         nsp += track->GetNcls(iDet);
1829       if (nsp) {
1830         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1831         track->SetTrackPointArray(sp);
1832         Int_t isptrack = 0;
1833         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1834           AliTracker *tracker = fTracker[iDet];
1835           if (!tracker) continue;
1836           Int_t nspdet = track->GetNcls(iDet);
1837           if (nspdet <= 0) continue;
1838           track->GetClusters(iDet,idx);
1839           AliTrackPoint p;
1840           Int_t isp = 0;
1841           Int_t isp2 = 0;
1842           while (isp < nspdet) {
1843             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1844             const Int_t kNTPCmax = 159;
1845             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
1846             if (!isvalid) continue;
1847             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1848           }
1849         }
1850       }
1851     }
1852   if (fTracker[3]){
1853     fTracker[3]->UnloadClusters();
1854     fLoader[3]->UnloadRecPoints();
1855   }
1856 }