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