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