STEER/AliReconstruction.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // class for running the reconstruction                                      //
  21 //                                                                           //
  22 // Clusters and tracks are created for all detectors and all events by       //
  23 // typing:                                                                   //
  24 //                                                                           //
  25 //   AliReconstruction rec;                                                  //
  26 //   rec.Run();                                                              //
  27 //                                                                           //
  28 // The Run method returns kTRUE in case of successful execution.             //
  29 //                                                                           //
  30 // If the input to the reconstruction are not simulated digits but raw data, //
  31 // this can be specified by an argument of the Run method or by the method   //
  32 //                                                                           //
  33 //   rec.SetInput("...");                                                    //
  34 //                                                                           //
  35 // The input formats and the corresponding argument are:                     //
  36 // - DDL raw data files: directory name, ends with "/"                       //
  37 // - raw data root file: root file name, extension ".root"                   //
  38 // - raw data DATE file: DATE file name, any other non-empty string          //
  39 // - MC root files     : empty string, default                               //
  40 //                                                                           //
  41 // By default all events are reconstructed. The reconstruction can be        //
  42 // limited to a range of events by giving the index of the first and the     //
  43 // last event as an argument to the Run method or by calling                 //
  44 //                                                                           //
  45 //   rec.SetEventRange(..., ...);                                            //
  46 //                                                                           //
  47 // The index -1 (default) can be used for the last event to indicate no      //
  48 // upper limit of the event range.                                           //
  49 //                                                                           //
  50 // The name of the galice file can be changed from the default               //
  51 // "galice.root" by passing it as argument to the AliReconstruction          //
  52 // constructor or by                                                         //
  53 //                                                                           //
  54 //   rec.SetGAliceFile("...");                                               //
  55 //                                                                           //
  56 // The local reconstruction can be switched on or off for individual         //
  57 // detectors by                                                              //
  58 //                                                                           //
  59 //   rec.SetRunLocalReconstruction("...");                                   //
  60 //                                                                           //
  61 // The argument is a (case sensitive) string with the names of the           //
  62 // detectors separated by a space. The special string "ALL" selects all      //
  63 // available detectors. This is the default.                                 //
  64 //                                                                           //
  65 // The reconstruction of the primary vertex position can be switched off by  //
  66 //                                                                           //
  67 //   rec.SetRunVertexFinder(kFALSE);                                         //
  68 //                                                                           //
  69 // The tracking and the creation of ESD tracks can be switched on for        //
  70 // selected detectors by                                                     //
  71 //                                                                           //
  72 //   rec.SetRunTracking("...");                                              //
  73 //                                                                           //
  74 // Uniform/nonuniform field tracking switches (default: uniform field)       //
  75 //                                                                           //
  76 //   rec.SetUniformFieldTracking();  ( rec.SetNonuniformFieldTracking(); )   //
  77 //                                                                           //
  78 // The filling of additional ESD information can be steered by               //
  79 //                                                                           //
  80 //   rec.SetFillESD("...");                                                  //
  81 //                                                                           //
  82 // Again, for both methods the string specifies the list of detectors.       //
  83 // The default is "ALL".                                                     //
  84 //                                                                           //
  85 // The call of the shortcut method                                           //
  86 //                                                                           //
  87 //   rec.SetRunReconstruction("...");                                        //
  88 //                                                                           //
  89 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
  90 // SetFillESD with the same detector selecting string as argument.           //
  91 //                                                                           //
  92 // The reconstruction requires digits or raw data as input. For the creation //
  93 // of digits and raw data have a look at the class AliSimulation.            //
  94 //                                                                           //
  95 // For debug purposes the method SetCheckPointLevel can be used. If the      //
  96 // argument is greater than 0, files with ESD events will be written after   //
  97 // selected steps of the reconstruction for each event:                      //
  98 //   level 1: after tracking and after filling of ESD (final)                //
  99 //   level 2: in addition after each tracking step                           //
 100 //   level 3: in addition after the filling of ESD for each detector         //
 101 // If a final check point file exists for an event, this event will be       //
 102 // skipped in the reconstruction. The tracking and the filling of ESD for    //
 103 // a detector will be skipped as well, if the corresponding check point      //
 104 // file exists. The ESD event will then be loaded from the file instead.     //
 105 //                                                                           //
 106 ///////////////////////////////////////////////////////////////////////////////
 107
 108 #include <TArrayF.h>
 109 #include <TFile.h>
 110 #include <TSystem.h>
 111 #include <TROOT.h>
 112 #include <TPluginManager.h>
 113 #include <TStopwatch.h>
 114 #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 "AliESDVertex.h"
 127 #include "AliTracker.h"
 128 #include "AliVertexer.h"
 129 #include "AliHeader.h"
 130 #include "AliGenEventHeader.h"
 131 #include "AliPID.h"
 132 #include "AliESDpid.h"
 133 #include "AliESDtrack.h"
 134
 135 #include "AliRunTag.h"
 136 //#include "AliLHCTag.h"
 137 #include "AliDetectorTag.h"
 138 #include "AliEventTag.h"
 139
 140 #include "AliTrackPointArray.h"
 141 #include "AliCDBManager.h"
 142
 143 ClassImp(AliReconstruction)
 144
 145
 146 //_____________________________________________________________________________
 147 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
 148
 149 //_____________________________________________________________________________
 150 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
 151                                      const char* name, const char* title) :
 152   TNamed(name, title),
 153
 154   fRunLocalReconstruction("ALL"),
 155   fUniformField(kTRUE),
 156   fRunVertexFinder(kTRUE),
 157   fRunHLTTracking(kFALSE),
 158   fRunTracking("ALL"),
 159   fFillESD("ALL"),
 160   fGAliceFileName(gAliceFilename),
 161   fInput(""),
 162   fFirstEvent(0),
 163   fLastEvent(-1),
 164   fStopOnError(kFALSE),
 165   fCheckPointLevel(0),
 166   fOptions(),
 167
 168   fRunLoader(NULL),
 169   fRawReader(NULL),
 170
 171   fVertexer(NULL),
 172
 173   fWriteAlignmentData(kFALSE),
 174   fCDBUri(cdbUri)
 175 {
 176 // create reconstruction object with default parameters
 177
 178   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 179     fReconstructor[iDet] = NULL;
 180     fLoader[iDet] = NULL;
 181     fTracker[iDet] = NULL;
 182   }
 183   AliPID pid;
 184   // Import TGeo geometry
 185   TString geom(gSystem->DirName(gAliceFilename));
 186   geom += "/geometry.root";
 187   TGeoManager::Import(geom.Data());
 188 }
 189
 190 //_____________________________________________________________________________
 191 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
 192   TNamed(rec),
 193
 194   fRunLocalReconstruction(rec.fRunLocalReconstruction),
 195   fUniformField(rec.fUniformField),
 196   fRunVertexFinder(rec.fRunVertexFinder),
 197   fRunHLTTracking(rec.fRunHLTTracking),
 198   fRunTracking(rec.fRunTracking),
 199   fFillESD(rec.fFillESD),
 200   fGAliceFileName(rec.fGAliceFileName),
 201   fInput(rec.fInput),
 202   fFirstEvent(rec.fFirstEvent),
 203   fLastEvent(rec.fLastEvent),
 204   fStopOnError(rec.fStopOnError),
 205   fCheckPointLevel(0),
 206   fOptions(),
 207
 208   fRunLoader(NULL),
 209   fRawReader(NULL),
 210
 211   fVertexer(NULL),
 212
 213   fWriteAlignmentData(rec.fWriteAlignmentData),
 214   fCDBUri(rec.fCDBUri)
 215 {
 216 // copy constructor
 217
 218   for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
 219     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
 220   }
 221   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 222     fReconstructor[iDet] = NULL;
 223     fLoader[iDet] = NULL;
 224     fTracker[iDet] = NULL;
 225   }
 226 }
 227
 228 //_____________________________________________________________________________
 229 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 230 {
 231 // assignment operator
 232
 233   this->~AliReconstruction();
 234   new(this) AliReconstruction(rec);
 235   return *this;
 236 }
 237
 238 //_____________________________________________________________________________
 239 AliReconstruction::~AliReconstruction()
 240 {
 241 // clean up
 242
 243   CleanUp();
 244   fOptions.Delete();
 245 }
 246
 247 //_____________________________________________________________________________
 248 void AliReconstruction::InitCDBStorage()
 249 {
 250 // activate a default CDB storage
 251 // First check if we have any CDB storage set, because it is used
 252 // to retrieve the calibration and alignment constants
 253
 254   AliCDBManager* man = AliCDBManager::Instance();
 255   if (!man->IsDefaultStorageSet())
 256   {
 257     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 258     AliWarning("Default CDB storage not yet set");
 259     AliWarning(Form("Using default storage declared in AliSimulation: %s",fCDBUri.Data()));
 260     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 261     SetDefaultStorage(fCDBUri);
 262   }
 263
 264 }
 265
 266 //_____________________________________________________________________________
 267 void AliReconstruction::SetDefaultStorage(const char* uri) {
 268 // activate a default CDB storage
 269
 270    AliCDBManager::Instance()->SetDefaultStorage(uri);
 271
 272 }
 273
 274 //_____________________________________________________________________________
 275 void AliReconstruction::SetSpecificStorage(const char* detName, const char* uri) {
 276 // activate a detector-specific CDB storage
 277
 278    AliCDBManager::Instance()->SetSpecificStorage(detName, uri);
 279
 280 }
 281
 282
 283 //_____________________________________________________________________________
 284 void AliReconstruction::SetGAliceFile(const char* fileName)
 285 {
 286 // set the name of the galice file
 287
 288   fGAliceFileName = fileName;
 289 }
 290
 291 //_____________________________________________________________________________
 292 void AliReconstruction::SetOption(const char* detector, const char* option)
 293 {
 294 // set options for the reconstruction of a detector
 295
 296   TObject* obj = fOptions.FindObject(detector);
 297   if (obj) fOptions.Remove(obj);
 298   fOptions.Add(new TNamed(detector, option));
 299 }
 300
 301
 302 //_____________________________________________________________________________
 303 Bool_t AliReconstruction::Run(const char* input,
 304                               Int_t firstEvent, Int_t lastEvent)
 305 {
 306 // run the reconstruction
 307
 308   InitCDBStorage();
 309
 310   // set the input
 311   if (!input) input = fInput.Data();
 312   TString fileName(input);
 313   if (fileName.EndsWith("/")) {
 314     fRawReader = new AliRawReaderFile(fileName);
 315   } else if (fileName.EndsWith(".root")) {
 316     fRawReader = new AliRawReaderRoot(fileName);
 317   } else if (!fileName.IsNull()) {
 318     fRawReader = new AliRawReaderDate(fileName);
 319     fRawReader->SelectEvents(7);
 320   }
 321
 322   // get the run loader
 323   if (!InitRunLoader()) return kFALSE;
 324
 325   // local reconstruction
 326   if (!fRunLocalReconstruction.IsNull()) {
 327     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
 328       if (fStopOnError) {CleanUp(); return kFALSE;}
 329     }
 330   }
 331 //  if (!fRunVertexFinder && fRunTracking.IsNull() &&
 332 //      fFillESD.IsNull()) return kTRUE;
 333
 334   // get vertexer
 335   if (fRunVertexFinder && !CreateVertexer()) {
 336     if (fStopOnError) {
 337       CleanUp();
 338       return kFALSE;
 339     }
 340   }
 341
 342   // get trackers
 343   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
 344     if (fStopOnError) {
 345       CleanUp();
 346       return kFALSE;
 347     }
 348   }
 349
 350
 351   TStopwatch stopwatch;
 352   stopwatch.Start();
 353
 354   // get the possibly already existing ESD file and tree
 355   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
 356   TFile* fileOld = NULL;
 357   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
 358   if (!gSystem->AccessPathName("AliESDs.root")){
 359     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
 360     fileOld = TFile::Open("AliESDs.old.root");
 361     if (fileOld && fileOld->IsOpen()) {
 362       treeOld = (TTree*) fileOld->Get("esdTree");
 363       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
 364       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
 365       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
 366     }
 367   }
 368
 369   // create the ESD output file and tree
 370   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 371   if (!file->IsOpen()) {
 372     AliError("opening AliESDs.root failed");
 373     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 374   }
 375   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 376   tree->Branch("ESD", "AliESD", &esd);
 377   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
 378   hlttree->Branch("ESD", "AliESD", &hltesd);
 379   delete esd; delete hltesd;
 380   esd = NULL; hltesd = NULL;
 381   gROOT->cd();
 382
 383   // loop over events
 384   if (fRawReader) fRawReader->RewindEvents();
 385
 386   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 387     if (fRawReader) fRawReader->NextEvent();
 388     if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
 389       // copy old ESD to the new one
 390       if (treeOld) {
 391         treeOld->SetBranchAddress("ESD", &esd);
 392         treeOld->GetEntry(iEvent);
 393       }
 394       tree->Fill();
 395       if (hlttreeOld) {
 396         hlttreeOld->SetBranchAddress("ESD", &hltesd);
 397         hlttreeOld->GetEntry(iEvent);
 398       }
 399       hlttree->Fill();
 400       continue;
 401     }
 402
 403     AliInfo(Form("processing event %d", iEvent));
 404     fRunLoader->GetEvent(iEvent);
 405
 406     char fileName[256];
 407     sprintf(fileName, "ESD_%d.%d_final.root",
 408             fRunLoader->GetHeader()->GetRun(),
 409             fRunLoader->GetHeader()->GetEventNrInRun());
 410     if (!gSystem->AccessPathName(fileName)) continue;
 411
 412     // local reconstruction
 413     if (!fRunLocalReconstruction.IsNull()) {
 414       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
 415         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 416       }
 417     }
 418
 419     esd = new AliESD; hltesd = new AliESD;
 420     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 421     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 422     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 423     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 424     if (gAlice) {
 425       esd->SetMagneticField(AliTracker::GetBz());
 426       hltesd->SetMagneticField(AliTracker::GetBz());
 427     } else {
 428       // ???
 429     }
 430
 431     // vertex finder
 432     if (fRunVertexFinder) {
 433       if (!ReadESD(esd, "vertex")) {
 434         if (!RunVertexFinder(esd)) {
 435           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 436         }
 437         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 438       }
 439     }
 440
 441     // HLT tracking
 442     if (!fRunTracking.IsNull()) {
 443       if (fRunHLTTracking) {
 444         hltesd->SetVertex(esd->GetVertex());
 445         if (!RunHLTTracking(hltesd)) {
 446           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 447         }
 448       }
 449     }
 450
 451     // barrel tracking
 452     if (!fRunTracking.IsNull()) {
 453       if (!ReadESD(esd, "tracking")) {
 454         if (!RunTracking(esd)) {
 455           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 456         }
 457         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 458       }
 459     }
 460
 461     // fill ESD
 462     if (!fFillESD.IsNull()) {
 463       if (!FillESD(esd, fFillESD)) {
 464         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 465       }
 466     }
 467
 468     // combined PID
 469     AliESDpid::MakePID(esd);
 470     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 471
 472     // write ESD
 473     tree->Fill();
 474     // write HLT ESD
 475     hlttree->Fill();
 476
 477     if (fCheckPointLevel > 0)  WriteESD(esd, "final");
 478
 479     delete esd; delete hltesd;
 480     esd = NULL; hltesd = NULL;
 481   }
 482
 483   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
 484                stopwatch.RealTime(),stopwatch.CpuTime()));
 485
 486   file->cd();
 487   tree->Write();
 488   hlttree->Write();
 489
 490   // Create tags for the events in the ESD tree (the ESD tree is always present)
 491   // In case of empty events the tags will contain dummy values
 492   CreateTag(file);
 493   CleanUp(file, fileOld);
 494
 495   return kTRUE;
 496 }
 497
 498
 499 //_____________________________________________________________________________
 500 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 501 {
 502 // run the local reconstruction
 503
 504   TStopwatch stopwatch;
 505   stopwatch.Start();
 506
 507   TString detStr = detectors;
 508   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 509     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 510     AliReconstructor* reconstructor = GetReconstructor(iDet);
 511     if (!reconstructor) continue;
 512     if (reconstructor->HasLocalReconstruction()) continue;
 513
 514     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 515     TStopwatch stopwatchDet;
 516     stopwatchDet.Start();
 517     if (fRawReader) {
 518       fRawReader->RewindEvents();
 519       reconstructor->Reconstruct(fRunLoader, fRawReader);
 520     } else {
 521       reconstructor->Reconstruct(fRunLoader);
 522     }
 523     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 524                  fgkDetectorName[iDet],
 525                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 526   }
 527
 528   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 529     AliError(Form("the following detectors were not found: %s",
 530                   detStr.Data()));
 531     if (fStopOnError) return kFALSE;
 532   }
 533
 534   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 535                stopwatch.RealTime(),stopwatch.CpuTime()));
 536
 537   return kTRUE;
 538 }
 539
 540 //_____________________________________________________________________________
 541 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 542 {
 543 // run the local reconstruction
 544
 545   TStopwatch stopwatch;
 546   stopwatch.Start();
 547
 548   TString detStr = detectors;
 549   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 550     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 551     AliReconstructor* reconstructor = GetReconstructor(iDet);
 552     if (!reconstructor) continue;
 553     AliLoader* loader = fLoader[iDet];
 554
 555     // conversion of digits
 556     if (fRawReader && reconstructor->HasDigitConversion()) {
 557       AliInfo(Form("converting raw data digits into root objects for %s",
 558                    fgkDetectorName[iDet]));
 559       TStopwatch stopwatchDet;
 560       stopwatchDet.Start();
 561       loader->LoadDigits("update");
 562       loader->CleanDigits();
 563       loader->MakeDigitsContainer();
 564       TTree* digitsTree = loader->TreeD();
 565       reconstructor->ConvertDigits(fRawReader, digitsTree);
 566       loader->WriteDigits("OVERWRITE");
 567       loader->UnloadDigits();
 568       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 569                    fgkDetectorName[iDet],
 570                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 571     }
 572
 573     // local reconstruction
 574     if (!reconstructor->HasLocalReconstruction()) continue;
 575     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 576     TStopwatch stopwatchDet;
 577     stopwatchDet.Start();
 578     loader->LoadRecPoints("update");
 579     loader->CleanRecPoints();
 580     loader->MakeRecPointsContainer();
 581     TTree* clustersTree = loader->TreeR();
 582     if (fRawReader && !reconstructor->HasDigitConversion()) {
 583       reconstructor->Reconstruct(fRawReader, clustersTree);
 584     } else {
 585       loader->LoadDigits("read");
 586       TTree* digitsTree = loader->TreeD();
 587       if (!digitsTree) {
 588         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 589         if (fStopOnError) return kFALSE;
 590       } else {
 591         reconstructor->Reconstruct(digitsTree, clustersTree);
 592       }
 593       loader->UnloadDigits();
 594     }
 595     loader->WriteRecPoints("OVERWRITE");
 596     loader->UnloadRecPoints();
 597     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
 598                     fgkDetectorName[iDet],
 599                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 600   }
 601
 602   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 603     AliError(Form("the following detectors were not found: %s",
 604                   detStr.Data()));
 605     if (fStopOnError) return kFALSE;
 606   }
 607
 608   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 609                stopwatch.RealTime(),stopwatch.CpuTime()));
 610
 611   return kTRUE;
 612 }
 613
 614 //_____________________________________________________________________________
 615 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 616 {
 617 // run the barrel tracking
 618
 619   TStopwatch stopwatch;
 620   stopwatch.Start();
 621
 622   AliESDVertex* vertex = NULL;
 623   Double_t vtxPos[3] = {0, 0, 0};
 624   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 625   TArrayF mcVertex(3);
 626   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 627     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 628     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 629   }
 630
 631   if (fVertexer) {
 632     AliInfo("running the ITS vertex finder");
 633     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 634     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 635     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 636     if(!vertex){
 637       AliWarning("Vertex not found");
 638       vertex = new AliESDVertex();
 639     }
 640     else {
 641       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 642     }
 643
 644   } else {
 645     AliInfo("getting the primary vertex from MC");
 646     vertex = new AliESDVertex(vtxPos, vtxErr);
 647   }
 648
 649   if (vertex) {
 650     vertex->GetXYZ(vtxPos);
 651     vertex->GetSigmaXYZ(vtxErr);
 652   } else {
 653     AliWarning("no vertex reconstructed");
 654     vertex = new AliESDVertex(vtxPos, vtxErr);
 655   }
 656   esd->SetVertex(vertex);
 657   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 658     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
 659   }
 660   delete vertex;
 661
 662   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 663                stopwatch.RealTime(),stopwatch.CpuTime()));
 664
 665   return kTRUE;
 666 }
 667
 668 //_____________________________________________________________________________
 669 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
 670 {
 671 // run the HLT barrel tracking
 672
 673   TStopwatch stopwatch;
 674   stopwatch.Start();
 675
 676   if (!fRunLoader) {
 677     AliError("Missing runLoader!");
 678     return kFALSE;
 679   }
 680
 681   AliInfo("running HLT tracking");
 682
 683   // Get a pointer to the HLT reconstructor
 684   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
 685   if (!reconstructor) return kFALSE;
 686
 687   // TPC + ITS
 688   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 689     TString detName = fgkDetectorName[iDet];
 690     AliDebug(1, Form("%s HLT tracking", detName.Data()));
 691     reconstructor->SetOption(detName.Data());
 692     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
 693     if (!tracker) {
 694       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
 695       if (fStopOnError) return kFALSE;
 696       continue;
 697     }
 698     Double_t vtxPos[3];
 699     Double_t vtxErr[3]={0.005,0.005,0.010};
 700     const AliESDVertex *vertex = esd->GetVertex();
 701     vertex->GetXYZ(vtxPos);
 702     tracker->SetVertex(vtxPos,vtxErr);
 703     if(iDet != 1) {
 704       fLoader[iDet]->LoadRecPoints("read");
 705       TTree* tree = fLoader[iDet]->TreeR();
 706       if (!tree) {
 707         AliError(Form("Can't get the %s cluster tree", detName.Data()));
 708         return kFALSE;
 709       }
 710       tracker->LoadClusters(tree);
 711     }
 712     if (tracker->Clusters2Tracks(esd) != 0) {
 713       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
 714       return kFALSE;
 715     }
 716     if(iDet != 1) {
 717       tracker->UnloadClusters();
 718     }
 719     delete tracker;
 720   }
 721
 722   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 723                stopwatch.RealTime(),stopwatch.CpuTime()));
 724
 725   return kTRUE;
 726 }
 727
 728 //_____________________________________________________________________________
 729 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
 730 {
 731 // run the barrel tracking
 732
 733   TStopwatch stopwatch;
 734   stopwatch.Start();
 735
 736   AliInfo("running tracking");
 737
 738   // pass 1: TPC + ITS inwards
 739   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 740     if (!fTracker[iDet]) continue;
 741     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
 742
 743     // load clusters
 744     fLoader[iDet]->LoadRecPoints("read");
 745     TTree* tree = fLoader[iDet]->TreeR();
 746     if (!tree) {
 747       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 748       return kFALSE;
 749     }
 750     fTracker[iDet]->LoadClusters(tree);
 751
 752     // run tracking
 753     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
 754       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
 755       return kFALSE;
 756     }
 757     if (fCheckPointLevel > 1) {
 758       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
 759     }
 760     // preliminary PID in TPC needed by the ITS tracker
 761     if (iDet == 1) {
 762       GetReconstructor(1)->FillESD(fRunLoader, esd);
 763       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
 764       AliESDpid::MakePID(esd);
 765     }
 766   }
 767
 768   // pass 2: ALL backwards
 769   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 770     if (!fTracker[iDet]) continue;
 771     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
 772
 773     // load clusters
 774     if (iDet > 1) {     // all except ITS, TPC
 775       TTree* tree = NULL;
 776       fLoader[iDet]->LoadRecPoints("read");
 777       tree = fLoader[iDet]->TreeR();
 778       if (!tree) {
 779         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 780         return kFALSE;
 781       }
 782       fTracker[iDet]->LoadClusters(tree);
 783     }
 784
 785     // run tracking
 786     if (fTracker[iDet]->PropagateBack(esd) != 0) {
 787       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
 788       return kFALSE;
 789     }
 790     if (fCheckPointLevel > 1) {
 791       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
 792     }
 793
 794     // unload clusters
 795     if (iDet > 2) {     // all except ITS, TPC, TRD
 796       fTracker[iDet]->UnloadClusters();
 797       fLoader[iDet]->UnloadRecPoints();
 798     }
 799     // updated PID in TPC needed by the ITS tracker -MI
 800     if (iDet == 1) {
 801       GetReconstructor(1)->FillESD(fRunLoader, esd);
 802       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
 803       AliESDpid::MakePID(esd);
 804     }
 805   }
 806
 807   // write space-points to the ESD in case alignment data output
 808   // is switched on
 809   if (fWriteAlignmentData)
 810     WriteAlignmentData(esd);
 811
 812   // pass 3: TRD + TPC + ITS refit inwards
 813   for (Int_t iDet = 2; iDet >= 0; iDet--) {
 814     if (!fTracker[iDet]) continue;
 815     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
 816
 817     // run tracking
 818     if (fTracker[iDet]->RefitInward(esd) != 0) {
 819       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
 820       return kFALSE;
 821     }
 822     if (fCheckPointLevel > 1) {
 823       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
 824     }
 825
 826     // unload clusters
 827     fTracker[iDet]->UnloadClusters();
 828     fLoader[iDet]->UnloadRecPoints();
 829   }
 830   //
 831   // Propagate track to the vertex - if not done by ITS
 832   //
 833   Int_t ntracks = esd->GetNumberOfTracks();
 834   for (Int_t itrack=0; itrack<ntracks; itrack++){
 835     const Double_t kRadius  = 3;   // beam pipe radius
 836     const Double_t kMaxStep = 5;   // max step
 837     const Double_t kMaxD    = 123456;  // max distance to prim vertex
 838     Double_t       fieldZ   = AliTracker::GetBz();  //
 839     AliESDtrack * track = esd->GetTrack(itrack);
 840     if (!track) continue;
 841     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
 842     track->PropagateTo(kRadius, track->GetMass(),kMaxStep,kTRUE);
 843     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
 844   }
 845
 846   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 847                stopwatch.RealTime(),stopwatch.CpuTime()));
 848
 849   return kTRUE;
 850 }
 851
 852 //_____________________________________________________________________________
 853 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
 854 {
 855 // fill the event summary data
 856
 857   TStopwatch stopwatch;
 858   stopwatch.Start();
 859   AliInfo("filling ESD");
 860
 861   TString detStr = detectors;
 862   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 863     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 864     AliReconstructor* reconstructor = GetReconstructor(iDet);
 865     if (!reconstructor) continue;
 866
 867     if (!ReadESD(esd, fgkDetectorName[iDet])) {
 868       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
 869       TTree* clustersTree = NULL;
 870       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 871         fLoader[iDet]->LoadRecPoints("read");
 872         clustersTree = fLoader[iDet]->TreeR();
 873         if (!clustersTree) {
 874           AliError(Form("Can't get the %s clusters tree",
 875                         fgkDetectorName[iDet]));
 876           if (fStopOnError) return kFALSE;
 877         }
 878       }
 879       if (fRawReader && !reconstructor->HasDigitConversion()) {
 880         reconstructor->FillESD(fRawReader, clustersTree, esd);
 881       } else {
 882         TTree* digitsTree = NULL;
 883         if (fLoader[iDet]) {
 884           fLoader[iDet]->LoadDigits("read");
 885           digitsTree = fLoader[iDet]->TreeD();
 886           if (!digitsTree) {
 887             AliError(Form("Can't get the %s digits tree",
 888                           fgkDetectorName[iDet]));
 889             if (fStopOnError) return kFALSE;
 890           }
 891         }
 892         reconstructor->FillESD(digitsTree, clustersTree, esd);
 893         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
 894       }
 895       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 896         fLoader[iDet]->UnloadRecPoints();
 897       }
 898
 899       if (fRawReader) {
 900         reconstructor->FillESD(fRunLoader, fRawReader, esd);
 901       } else {
 902         reconstructor->FillESD(fRunLoader, esd);
 903       }
 904       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
 905     }
 906   }
 907
 908   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 909     AliError(Form("the following detectors were not found: %s",
 910                   detStr.Data()));
 911     if (fStopOnError) return kFALSE;
 912   }
 913
 914   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 915                stopwatch.RealTime(),stopwatch.CpuTime()));
 916
 917   return kTRUE;
 918 }
 919
 920
 921 //_____________________________________________________________________________
 922 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
 923 {
 924 // check whether detName is contained in detectors
 925 // if yes, it is removed from detectors
 926
 927   // check if all detectors are selected
 928   if ((detectors.CompareTo("ALL") == 0) ||
 929       detectors.BeginsWith("ALL ") ||
 930       detectors.EndsWith(" ALL") ||
 931       detectors.Contains(" ALL ")) {
 932     detectors = "ALL";
 933     return kTRUE;
 934   }
 935
 936   // search for the given detector
 937   Bool_t result = kFALSE;
 938   if ((detectors.CompareTo(detName) == 0) ||
 939       detectors.BeginsWith(detName+" ") ||
 940       detectors.EndsWith(" "+detName) ||
 941       detectors.Contains(" "+detName+" ")) {
 942     detectors.ReplaceAll(detName, "");
 943     result = kTRUE;
 944   }
 945
 946   // clean up the detectors string
 947   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
 948   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
 949   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
 950
 951   return result;
 952 }
 953
 954 //_____________________________________________________________________________
 955 Bool_t AliReconstruction::InitRunLoader()
 956 {
 957 // get or create the run loader
 958
 959   if (gAlice) delete gAlice;
 960   gAlice = NULL;
 961
 962   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
 963     // load all base libraries to get the loader classes
 964     TString libs = gSystem->GetLibraries();
 965     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 966       TString detName = fgkDetectorName[iDet];
 967       if (detName == "HLT") continue;
 968       if (libs.Contains("lib" + detName + "base.so")) continue;
 969       gSystem->Load("lib" + detName + "base.so");
 970     }
 971     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
 972     if (!fRunLoader) {
 973       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
 974       CleanUp();
 975       return kFALSE;
 976     }
 977     fRunLoader->CdGAFile();
 978     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
 979       if (fRunLoader->LoadgAlice() == 0) {
 980         gAlice = fRunLoader->GetAliRun();
 981         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
 982       }
 983     }
 984     if (!gAlice && !fRawReader) {
 985       AliError(Form("no gAlice object found in file %s",
 986                     fGAliceFileName.Data()));
 987       CleanUp();
 988       return kFALSE;
 989     }
 990
 991   } else {               // galice.root does not exist
 992     if (!fRawReader) {
 993       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
 994       CleanUp();
 995       return kFALSE;
 996     }
 997     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
 998                                     AliConfig::GetDefaultEventFolderName(),
 999                                     "recreate");
1000     if (!fRunLoader) {
1001       AliError(Form("could not create run loader in file %s",
1002                     fGAliceFileName.Data()));
1003       CleanUp();
1004       return kFALSE;
1005     }
1006     fRunLoader->MakeTree("E");
1007     Int_t iEvent = 0;
1008     while (fRawReader->NextEvent()) {
1009       fRunLoader->SetEventNumber(iEvent);
1010       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1011                                      iEvent, iEvent);
1012       fRunLoader->MakeTree("H");
1013       fRunLoader->TreeE()->Fill();
1014       iEvent++;
1015     }
1016     fRawReader->RewindEvents();
1017     fRunLoader->WriteHeader("OVERWRITE");
1018     fRunLoader->CdGAFile();
1019     fRunLoader->Write(0, TObject::kOverwrite);
1020 //    AliTracker::SetFieldMap(???);
1021   }
1022
1023   return kTRUE;
1024 }
1025
1026 //_____________________________________________________________________________
1027 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1028 {
1029 // get the reconstructor object and the loader for a detector
1030
1031   if (fReconstructor[iDet]) return fReconstructor[iDet];
1032
1033   // load the reconstructor object
1034   TPluginManager* pluginManager = gROOT->GetPluginManager();
1035   TString detName = fgkDetectorName[iDet];
1036   TString recName = "Ali" + detName + "Reconstructor";
1037   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1038
1039   if (detName == "HLT") {
1040     if (!gROOT->GetClass("AliLevel3")) {
1041       gSystem->Load("libAliL3Src.so");
1042       gSystem->Load("libAliL3Misc.so");
1043       gSystem->Load("libAliL3Hough.so");
1044       gSystem->Load("libAliL3Comp.so");
1045     }
1046   }
1047
1048   AliReconstructor* reconstructor = NULL;
1049   // first check if a plugin is defined for the reconstructor
1050   TPluginHandler* pluginHandler =
1051     pluginManager->FindHandler("AliReconstructor", detName);
1052   // if not, add a plugin for it
1053   if (!pluginHandler) {
1054     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1055     TString libs = gSystem->GetLibraries();
1056     if (libs.Contains("lib" + detName + "base.so") ||
1057         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1058       pluginManager->AddHandler("AliReconstructor", detName,
1059                                 recName, detName + "rec", recName + "()");
1060     } else {
1061       pluginManager->AddHandler("AliReconstructor", detName,
1062                                 recName, detName, recName + "()");
1063     }
1064     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1065   }
1066   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1067     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1068   }
1069   if (reconstructor) {
1070     TObject* obj = fOptions.FindObject(detName.Data());
1071     if (obj) reconstructor->SetOption(obj->GetTitle());
1072     reconstructor->Init(fRunLoader);
1073     fReconstructor[iDet] = reconstructor;
1074   }
1075
1076   // get or create the loader
1077   if (detName != "HLT") {
1078     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1079     if (!fLoader[iDet]) {
1080       AliConfig::Instance()
1081         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1082                                 detName, detName);
1083       // first check if a plugin is defined for the loader
1084       TPluginHandler* pluginHandler =
1085         pluginManager->FindHandler("AliLoader", detName);
1086       // if not, add a plugin for it
1087       if (!pluginHandler) {
1088         TString loaderName = "Ali" + detName + "Loader";
1089         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1090         pluginManager->AddHandler("AliLoader", detName,
1091                                   loaderName, detName + "base",
1092                                   loaderName + "(const char*, TFolder*)");
1093         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1094       }
1095       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1096         fLoader[iDet] =
1097           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1098                                                  fRunLoader->GetEventFolder());
1099       }
1100       if (!fLoader[iDet]) {   // use default loader
1101         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1102       }
1103       if (!fLoader[iDet]) {
1104         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1105         if (fStopOnError) return NULL;
1106       } else {
1107         fRunLoader->AddLoader(fLoader[iDet]);
1108         fRunLoader->CdGAFile();
1109         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1110         fRunLoader->Write(0, TObject::kOverwrite);
1111       }
1112     }
1113   }
1114
1115   return reconstructor;
1116 }
1117
1118 //_____________________________________________________________________________
1119 Bool_t AliReconstruction::CreateVertexer()
1120 {
1121 // create the vertexer
1122
1123   fVertexer = NULL;
1124   AliReconstructor* itsReconstructor = GetReconstructor(0);
1125   if (itsReconstructor) {
1126     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1127   }
1128   if (!fVertexer) {
1129     AliWarning("couldn't create a vertexer for ITS");
1130     if (fStopOnError) return kFALSE;
1131   }
1132
1133   return kTRUE;
1134 }
1135
1136 //_____________________________________________________________________________
1137 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1138 {
1139 // create the trackers
1140
1141   TString detStr = detectors;
1142   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1143     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1144     AliReconstructor* reconstructor = GetReconstructor(iDet);
1145     if (!reconstructor) continue;
1146     TString detName = fgkDetectorName[iDet];
1147     if (detName == "HLT") {
1148       fRunHLTTracking = kTRUE;
1149       continue;
1150     }
1151
1152     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1153     if (!fTracker[iDet] && (iDet < 7)) {
1154       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1155       if (fStopOnError) return kFALSE;
1156     }
1157   }
1158
1159   return kTRUE;
1160 }
1161
1162 //_____________________________________________________________________________
1163 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1164 {
1165 // delete trackers and the run loader and close and delete the file
1166
1167   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1168     delete fReconstructor[iDet];
1169     fReconstructor[iDet] = NULL;
1170     fLoader[iDet] = NULL;
1171     delete fTracker[iDet];
1172     fTracker[iDet] = NULL;
1173   }
1174   delete fVertexer;
1175   fVertexer = NULL;
1176
1177   delete fRunLoader;
1178   fRunLoader = NULL;
1179   delete fRawReader;
1180   fRawReader = NULL;
1181
1182   if (file) {
1183     file->Close();
1184     delete file;
1185   }
1186
1187   if (fileOld) {
1188     fileOld->Close();
1189     delete fileOld;
1190     gSystem->Unlink("AliESDs.old.root");
1191   }
1192 }
1193
1194
1195 //_____________________________________________________________________________
1196 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1197 {
1198 // read the ESD event from a file
1199
1200   if (!esd) return kFALSE;
1201   char fileName[256];
1202   sprintf(fileName, "ESD_%d.%d_%s.root",
1203           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1204   if (gSystem->AccessPathName(fileName)) return kFALSE;
1205
1206   AliInfo(Form("reading ESD from file %s", fileName));
1207   AliDebug(1, Form("reading ESD from file %s", fileName));
1208   TFile* file = TFile::Open(fileName);
1209   if (!file || !file->IsOpen()) {
1210     AliError(Form("opening %s failed", fileName));
1211     delete file;
1212     return kFALSE;
1213   }
1214
1215   gROOT->cd();
1216   delete esd;
1217   esd = (AliESD*) file->Get("ESD");
1218   file->Close();
1219   delete file;
1220   return kTRUE;
1221 }
1222
1223 //_____________________________________________________________________________
1224 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1225 {
1226 // write the ESD event to a file
1227
1228   if (!esd) return;
1229   char fileName[256];
1230   sprintf(fileName, "ESD_%d.%d_%s.root",
1231           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1232
1233   AliDebug(1, Form("writing ESD to file %s", fileName));
1234   TFile* file = TFile::Open(fileName, "recreate");
1235   if (!file || !file->IsOpen()) {
1236     AliError(Form("opening %s failed", fileName));
1237   } else {
1238     esd->Write("ESD");
1239     file->Close();
1240   }
1241   delete file;
1242 }
1243
1244
1245
1246
1247 //_____________________________________________________________________________
1248 void AliReconstruction::CreateTag(TFile* file)
1249 {
1250   /////////////
1251   //muon code//
1252   ////////////
1253   Double_t fMUONMASS = 0.105658369;
1254   //Variables
1255   Double_t fX,fY,fZ ;
1256   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1257   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1258   Int_t fCharge;
1259   TLorentzVector fEPvector;
1260
1261   Float_t fZVertexCut = 10.0;
1262   Float_t fRhoVertexCut = 2.0;
1263
1264   Float_t fLowPtCut = 1.0;
1265   Float_t fHighPtCut = 3.0;
1266   Float_t fVeryHighPtCut = 10.0;
1267   ////////////
1268
1269   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1270
1271   // Creates the tags for all the events in a given ESD file
1272   Int_t ntrack;
1273   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1274   Int_t nPos, nNeg, nNeutr;
1275   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1276   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1277   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1278   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1279   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1280
1281   AliRunTag *tag = new AliRunTag();
1282   AliEventTag *evTag = new AliEventTag();
1283   TTree ttag("T","A Tree with event tags");
1284   TBranch * btag = ttag.Branch("AliTAG", &tag);
1285   btag->SetCompressionLevel(9);
1286
1287   AliInfo(Form("Creating the tags......."));
1288
1289   if (!file || !file->IsOpen()) {
1290     AliError(Form("opening failed"));
1291     delete file;
1292     return ;
1293   }
1294   Int_t firstEvent = 0,lastEvent = 0;
1295   TTree *t = (TTree*) file->Get("esdTree");
1296   TBranch * b = t->GetBranch("ESD");
1297   AliESD *esd = 0;
1298   b->SetAddress(&esd);
1299
1300   tag->SetRunId(esd->GetRunNumber());
1301
1302   Int_t iNumberOfEvents = b->GetEntries();
1303   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1304     ntrack = 0;
1305     nPos = 0;
1306     nNeg = 0;
1307     nNeutr =0;
1308     nK0s = 0;
1309     nNeutrons = 0;
1310     nPi0s = 0;
1311     nGamas = 0;
1312     nProtons = 0;
1313     nKaons = 0;
1314     nPions = 0;
1315     nMuons = 0;
1316     nElectrons = 0;
1317     nCh1GeV = 0;
1318     nCh3GeV = 0;
1319     nCh10GeV = 0;
1320     nMu1GeV = 0;
1321     nMu3GeV = 0;
1322     nMu10GeV = 0;
1323     nEl1GeV = 0;
1324     nEl3GeV = 0;
1325     nEl10GeV = 0;
1326     maxPt = .0;
1327     meanPt = .0;
1328     totalP = .0;
1329
1330     b->GetEntry(iEventNumber);
1331     const AliESDVertex * vertexIn = esd->GetVertex();
1332
1333     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1334       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1335       UInt_t status = esdTrack->GetStatus();
1336
1337       //select only tracks with ITS refit
1338       if ((status&AliESDtrack::kITSrefit)==0) continue;
1339       //select only tracks with TPC refit
1340       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1341
1342       //select only tracks with the "combined PID"
1343       if ((status&AliESDtrack::kESDpid)==0) continue;
1344       Double_t p[3];
1345       esdTrack->GetPxPyPz(p);
1346       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1347       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1348       totalP += momentum;
1349       meanPt += fPt;
1350       if(fPt > maxPt) maxPt = fPt;
1351
1352       if(esdTrack->GetSign() > 0) {
1353         nPos++;
1354         if(fPt > fLowPtCut) nCh1GeV++;
1355         if(fPt > fHighPtCut) nCh3GeV++;
1356         if(fPt > fVeryHighPtCut) nCh10GeV++;
1357       }
1358       if(esdTrack->GetSign() < 0) {
1359         nNeg++;
1360         if(fPt > fLowPtCut) nCh1GeV++;
1361         if(fPt > fHighPtCut) nCh3GeV++;
1362         if(fPt > fVeryHighPtCut) nCh10GeV++;
1363       }
1364       if(esdTrack->GetSign() == 0) nNeutr++;
1365
1366       //PID
1367       Double_t prob[5];
1368       esdTrack->GetESDpid(prob);
1369
1370       Double_t rcc = 0.0;
1371       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1372       if(rcc == 0.0) continue;
1373       //Bayes' formula
1374       Double_t w[5];
1375       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1376
1377       //protons
1378       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1379       //kaons
1380       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1381       //pions
1382       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1383       //electrons
1384       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1385         nElectrons++;
1386         if(fPt > fLowPtCut) nEl1GeV++;
1387         if(fPt > fHighPtCut) nEl3GeV++;
1388         if(fPt > fVeryHighPtCut) nEl10GeV++;
1389       }
1390       ntrack++;
1391     }//track loop
1392
1393     /////////////
1394     //muon code//
1395     ////////////
1396     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1397     // loop over all reconstructed tracks (also first track of combination)
1398     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1399       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1400       if (muonTrack == 0x0) continue;
1401
1402       // Coordinates at vertex
1403       fZ = muonTrack->GetZ();
1404       fY = muonTrack->GetBendingCoor();
1405       fX = muonTrack->GetNonBendingCoor();
1406
1407       fThetaX = muonTrack->GetThetaX();
1408       fThetaY = muonTrack->GetThetaY();
1409
1410       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1411       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1412       fPxRec = fPzRec * TMath::Tan(fThetaX);
1413       fPyRec = fPzRec * TMath::Tan(fThetaY);
1414       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1415
1416       //ChiSquare of the track if needed
1417       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1418       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1419       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1420
1421       // total number of muons inside a vertex cut
1422       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1423         nMuons++;
1424         if(fEPvector.Pt() > fLowPtCut) {
1425           nMu1GeV++;
1426           if(fEPvector.Pt() > fHighPtCut) {
1427             nMu3GeV++;
1428             if (fEPvector.Pt() > fVeryHighPtCut) {
1429               nMu10GeV++;
1430             }
1431           }
1432         }
1433       }
1434     }//muon track loop
1435
1436     // Fill the event tags
1437     if(ntrack != 0)
1438       meanPt = meanPt/ntrack;
1439
1440     evTag->SetEventId(iEventNumber+1);
1441     evTag->SetVertexX(vertexIn->GetXv());
1442     evTag->SetVertexY(vertexIn->GetYv());
1443     evTag->SetVertexZ(vertexIn->GetZv());
1444
1445     evTag->SetT0VertexZ(esd->GetT0zVertex());
1446
1447     evTag->SetTrigger(esd->GetTrigger());
1448
1449     evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
1450     evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
1451     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1452     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1453
1454
1455     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1456     evTag->SetNumOfPosTracks(nPos);
1457     evTag->SetNumOfNegTracks(nNeg);
1458     evTag->SetNumOfNeutrTracks(nNeutr);
1459
1460     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1461     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1462     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1463     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1464
1465     evTag->SetNumOfProtons(nProtons);
1466     evTag->SetNumOfKaons(nKaons);
1467     evTag->SetNumOfPions(nPions);
1468     evTag->SetNumOfMuons(nMuons);
1469     evTag->SetNumOfElectrons(nElectrons);
1470     evTag->SetNumOfPhotons(nGamas);
1471     evTag->SetNumOfPi0s(nPi0s);
1472     evTag->SetNumOfNeutrons(nNeutrons);
1473     evTag->SetNumOfKaon0s(nK0s);
1474
1475     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1476     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1477     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1478     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1479     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1480     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1481     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1482     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1483     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1484
1485     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1486     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1487
1488     evTag->SetTotalMomentum(totalP);
1489     evTag->SetMeanPt(meanPt);
1490     evTag->SetMaxPt(maxPt);
1491
1492     tag->AddEventTag(*evTag);
1493   }
1494   lastEvent = iNumberOfEvents;
1495
1496   ttag.Fill();
1497   tag->Clear();
1498
1499   char fileName[256];
1500   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1501           tag->GetRunId(),firstEvent,lastEvent );
1502   AliInfo(Form("writing tags to file %s", fileName));
1503   AliDebug(1, Form("writing tags to file %s", fileName));
1504
1505   TFile* ftag = TFile::Open(fileName, "recreate");
1506   ftag->cd();
1507   ttag.Write();
1508   ftag->Close();
1509   file->cd();
1510   delete tag;
1511   delete evTag;
1512 }
1513
1514 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1515 {
1516   // Write space-points which are then used in the alignment procedures
1517   // For the moment only ITS, TRD and TPC
1518
1519   // Load TOF clusters
1520   if (fTracker[3]){
1521     fLoader[3]->LoadRecPoints("read");
1522     TTree* tree = fLoader[3]->TreeR();
1523     if (!tree) {
1524       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1525       return;
1526     }
1527     fTracker[3]->LoadClusters(tree);
1528   }
1529   Int_t ntracks = esd->GetNumberOfTracks();
1530   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1531     {
1532       AliESDtrack *track = esd->GetTrack(itrack);
1533       Int_t nsp = 0;
1534       UInt_t idx[200];
1535       for (Int_t iDet = 3; iDet >= 0; iDet--)
1536         nsp += track->GetNcls(iDet);
1537       if (nsp) {
1538         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1539         track->SetTrackPointArray(sp);
1540         Int_t isptrack = 0;
1541         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1542           AliTracker *tracker = fTracker[iDet];
1543           if (!tracker) continue;
1544           Int_t nspdet = track->GetNcls(iDet);
1545           if (nspdet <= 0) continue;
1546           track->GetClusters(iDet,idx);
1547           AliTrackPoint p;
1548           Int_t isp = 0;
1549           Int_t isp2 = 0;
1550           while (isp < nspdet) {
1551             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1552             if (!isvalid) continue;
1553             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1554           }
1555         }
1556       }
1557     }
1558   if (fTracker[3]){
1559     fTracker[3]->UnloadClusters();
1560     fLoader[3]->UnloadRecPoints();
1561   }
1562 }