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 // In case of raw-data reconstruction the user can modify the default        //
  51 // number of events per digits/clusters/tracks file. In case the option      //
  52 // is not used the number is set 1. In case the user provides 0, than        //
  53 // the number of events is equal to the number of events inside the          //
  54 // raw-data file (i.e. one digits/clusters/tracks file):                     //
  55 //                                                                           //
  56 //   rec.SetNumberOfEventsPerFile(...);                                      //
  57 //                                                                           //
  58 //                                                                           //
  59 // The name of the galice file can be changed from the default               //
  60 // "galice.root" by passing it as argument to the AliReconstruction          //
  61 // constructor or by                                                         //
  62 //                                                                           //
  63 //   rec.SetGAliceFile("...");                                               //
  64 //                                                                           //
  65 // The local reconstruction can be switched on or off for individual         //
  66 // detectors by                                                              //
  67 //                                                                           //
  68 //   rec.SetRunLocalReconstruction("...");                                   //
  69 //                                                                           //
  70 // The argument is a (case sensitive) string with the names of the           //
  71 // detectors separated by a space. The special string "ALL" selects all      //
  72 // available detectors. This is the default.                                 //
  73 //                                                                           //
  74 // The reconstruction of the primary vertex position can be switched off by  //
  75 //                                                                           //
  76 //   rec.SetRunVertexFinder(kFALSE);                                         //
  77 //                                                                           //
  78 // The tracking and the creation of ESD tracks can be switched on for        //
  79 // selected detectors by                                                     //
  80 //                                                                           //
  81 //   rec.SetRunTracking("...");                                              //
  82 //                                                                           //
  83 // Uniform/nonuniform field tracking switches (default: uniform field)       //
  84 //                                                                           //
  85 //   rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
  86 //                                                                           //
  87 // The filling of additional ESD information can be steered by               //
  88 //                                                                           //
  89 //   rec.SetFillESD("...");                                                  //
  90 //                                                                           //
  91 // Again, for both methods the string specifies the list of detectors.       //
  92 // The default is "ALL".                                                     //
  93 //                                                                           //
  94 // The call of the shortcut method                                           //
  95 //                                                                           //
  96 //   rec.SetRunReconstruction("...");                                        //
  97 //                                                                           //
  98 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
  99 // SetFillESD with the same detector selecting string as argument.           //
 100 //                                                                           //
 101 // The reconstruction requires digits or raw data as input. For the creation //
 102 // of digits and raw data have a look at the class AliSimulation.            //
 103 //                                                                           //
 104 // For debug purposes the method SetCheckPointLevel can be used. If the      //
 105 // argument is greater than 0, files with ESD events will be written after   //
 106 // selected steps of the reconstruction for each event:                      //
 107 //   level 1: after tracking and after filling of ESD (final)                //
 108 //   level 2: in addition after each tracking step                           //
 109 //   level 3: in addition after the filling of ESD for each detector         //
 110 // If a final check point file exists for an event, this event will be       //
 111 // skipped in the reconstruction. The tracking and the filling of ESD for    //
 112 // a detector will be skipped as well, if the corresponding check point      //
 113 // file exists. The ESD event will then be loaded from the file instead.     //
 114 //                                                                           //
 115 ///////////////////////////////////////////////////////////////////////////////
 116
 117 #include <TArrayF.h>
 118 #include <TFile.h>
 119 #include <TSystem.h>
 120 #include <TROOT.h>
 121 #include <TPluginManager.h>
 122 #include <TGeoManager.h>
 123 #include <TLorentzVector.h>
 124
 125 #include "AliReconstruction.h"
 126 #include "AliCodeTimer.h"
 127 #include "AliReconstructor.h"
 128 #include "AliLog.h"
 129 #include "AliRunLoader.h"
 130 #include "AliRun.h"
 131 #include "AliRawReaderFile.h"
 132 #include "AliRawReaderDate.h"
 133 #include "AliRawReaderRoot.h"
 134 #include "AliRawEventHeaderBase.h"
 135 #include "AliESDEvent.h"
 136 #include "AliESDMuonTrack.h"
 137 #include "AliESDfriend.h"
 138 #include "AliESDVertex.h"
 139 #include "AliESDcascade.h"
 140 #include "AliESDkink.h"
 141 #include "AliESDtrack.h"
 142 #include "AliESDCaloCluster.h"
 143 #include "AliMultiplicity.h"
 144 #include "AliTracker.h"
 145 #include "AliVertexer.h"
 146 #include "AliVertexerTracks.h"
 147 #include "AliV0vertexer.h"
 148 #include "AliCascadeVertexer.h"
 149 #include "AliHeader.h"
 150 #include "AliGenEventHeader.h"
 151 #include "AliPID.h"
 152 #include "AliESDpid.h"
 153 #include "AliESDtrack.h"
 154
 155 #include "AliESDTagCreator.h"
 156 #include "AliAODTagCreator.h"
 157
 158 #include "AliGeomManager.h"
 159 #include "AliTrackPointArray.h"
 160 #include "AliCDBManager.h"
 161 #include "AliCDBEntry.h"
 162 #include "AliAlignObj.h"
 163
 164 #include "AliCentralTrigger.h"
 165 #include "AliCTPRawStream.h"
 166
 167 #include "AliAODEvent.h"
 168 #include "AliAODHeader.h"
 169 #include "AliAODTrack.h"
 170 #include "AliAODVertex.h"
 171 #include "AliAODCluster.h"
 172
 173 #include "AliQualAssDataMaker.h"
 174
 175 ClassImp(AliReconstruction)
 176
 177
 178 //_____________________________________________________________________________
 179 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
 180
 181 //_____________________________________________________________________________
 182 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
 183                                      const char* name, const char* title) :
 184   TNamed(name, title),
 185
 186   fUniformField(kTRUE),
 187   fRunVertexFinder(kTRUE),
 188   fRunHLTTracking(kFALSE),
 189   fRunMuonTracking(kFALSE),
 190   fRunV0Finder(kTRUE),
 191   fRunCascadeFinder(kTRUE),
 192   fStopOnError(kFALSE),
 193   fWriteAlignmentData(kFALSE),
 194   fCleanESD(kTRUE),
 195   fWriteESDfriend(kFALSE),
 196   fWriteAOD(kFALSE),
 197   fFillTriggerESD(kTRUE),
 198
 199   fRunLocalReconstruction("ALL"),
 200   fRunTracking("ALL"),
 201   fFillESD("ALL"),
 202   fGAliceFileName(gAliceFilename),
 203   fInput(""),
 204   fEquipIdMap(""),
 205   fFirstEvent(0),
 206   fLastEvent(-1),
 207   fNumberOfEventsPerFile(1),
 208   fCheckPointLevel(0),
 209   fOptions(),
 210   fLoadAlignFromCDB(kTRUE),
 211   fLoadAlignData("ALL"),
 212   fESDPar(""),
 213
 214   fRunLoader(NULL),
 215   fRawReader(NULL),
 216
 217   fVertexer(NULL),
 218   fDiamondProfile(NULL),
 219
 220   fAlignObjArray(NULL),
 221   fCDBUri(cdbUri),
 222   fRemoteCDBUri(""),
 223   fSpecCDBUri()
 224 {
 225 // create reconstruction object with default parameters
 226
 227   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 228     fReconstructor[iDet] = NULL;
 229     fLoader[iDet] = NULL;
 230     fTracker[iDet] = NULL;
 231     fQualAssDataMaker[iDet] = NULL;
 232   }
 233   AliPID pid;
 234 }
 235
 236 //_____________________________________________________________________________
 237 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
 238   TNamed(rec),
 239
 240   fUniformField(rec.fUniformField),
 241   fRunVertexFinder(rec.fRunVertexFinder),
 242   fRunHLTTracking(rec.fRunHLTTracking),
 243   fRunMuonTracking(rec.fRunMuonTracking),
 244   fRunV0Finder(rec.fRunV0Finder),
 245   fRunCascadeFinder(rec.fRunCascadeFinder),
 246   fStopOnError(rec.fStopOnError),
 247   fWriteAlignmentData(rec.fWriteAlignmentData),
 248   fCleanESD(rec.fCleanESD),
 249   fWriteESDfriend(rec.fWriteESDfriend),
 250   fWriteAOD(rec.fWriteAOD),
 251   fFillTriggerESD(rec.fFillTriggerESD),
 252
 253   fRunLocalReconstruction(rec.fRunLocalReconstruction),
 254   fRunTracking(rec.fRunTracking),
 255   fFillESD(rec.fFillESD),
 256   fGAliceFileName(rec.fGAliceFileName),
 257   fInput(rec.fInput),
 258   fEquipIdMap(rec.fEquipIdMap),
 259   fFirstEvent(rec.fFirstEvent),
 260   fLastEvent(rec.fLastEvent),
 261   fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
 262   fCheckPointLevel(0),
 263   fOptions(),
 264   fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
 265   fLoadAlignData(rec.fLoadAlignData),
 266   fESDPar(rec.fESDPar),
 267
 268   fRunLoader(NULL),
 269   fRawReader(NULL),
 270
 271   fVertexer(NULL),
 272   fDiamondProfile(NULL),
 273
 274   fAlignObjArray(rec.fAlignObjArray),
 275   fCDBUri(rec.fCDBUri),
 276   fRemoteCDBUri(rec.fRemoteCDBUri),
 277   fSpecCDBUri()
 278 {
 279 // copy constructor
 280
 281   for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
 282     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
 283   }
 284   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 285     fReconstructor[iDet] = NULL;
 286     fLoader[iDet] = NULL;
 287     fTracker[iDet] = NULL;
 288     fQualAssDataMaker[iDet] = NULL;
 289   }
 290   for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
 291     if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
 292   }
 293 }
 294
 295 //_____________________________________________________________________________
 296 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 297 {
 298 // assignment operator
 299
 300   this->~AliReconstruction();
 301   new(this) AliReconstruction(rec);
 302   return *this;
 303 }
 304
 305 //_____________________________________________________________________________
 306 AliReconstruction::~AliReconstruction()
 307 {
 308 // clean up
 309
 310   CleanUp();
 311   fOptions.Delete();
 312   fSpecCDBUri.Delete();
 313
 314   AliCodeTimer::Instance()->Print();
 315 }
 316
 317 //_____________________________________________________________________________
 318 void AliReconstruction::InitCDBStorage()
 319 {
 320 // activate a default CDB storage
 321 // First check if we have any CDB storage set, because it is used
 322 // to retrieve the calibration and alignment constants
 323
 324   AliCDBManager* man = AliCDBManager::Instance();
 325   if (man->IsDefaultStorageSet())
 326   {
 327     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 328     AliWarning("Default CDB storage has been already set !");
 329     AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
 330     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 331     fCDBUri = "";
 332   }
 333   else {
 334     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 335     AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
 336     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 337     man->SetDefaultStorage(fCDBUri);
 338   }
 339
 340   // Remote storage (the Grid storage) is used if it is activated
 341   // and if the object is not found in the default storage
 342   if (man->IsRemoteStorageSet())
 343   {
 344     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 345     AliWarning("Remote CDB storage has been already set !");
 346     AliWarning(Form("Ignoring the remote storage declared in AliReconstruction: %s",fRemoteCDBUri.Data()));
 347     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 348     fRemoteCDBUri = "";
 349   }
 350   else {
 351     AliDebug(2,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 352     AliDebug(2, Form("Remote CDB storage is set to: %s",fRemoteCDBUri.Data()));
 353     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 354     man->SetRemoteStorage(fRemoteCDBUri);
 355   }
 356
 357   // Now activate the detector specific CDB storage locations
 358   for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
 359     TObject* obj = fSpecCDBUri[i];
 360     if (!obj) continue;
 361     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 362     AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
 363     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 364     man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
 365   }
 366   man->Print();
 367 }
 368
 369 //_____________________________________________________________________________
 370 void AliReconstruction::SetDefaultStorage(const char* uri) {
 371 // Store the desired default CDB storage location
 372 // Activate it later within the Run() method
 373
 374   fCDBUri = uri;
 375
 376 }
 377
 378 //_____________________________________________________________________________
 379 void AliReconstruction::SetRemoteStorage(const char* uri) {
 380 // Store the desired remote CDB storage location
 381 // Activate it later within the Run() method
 382 // Remote storage (the Grid storage) is used if it is activated
 383 // and if the object is not found in the default storage
 384
 385   fRemoteCDBUri = uri;
 386
 387 }
 388
 389 //_____________________________________________________________________________
 390 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
 391 // Store a detector-specific CDB storage location
 392 // Activate it later within the Run() method
 393
 394   AliCDBPath aPath(calibType);
 395   if(!aPath.IsValid()){
 396         // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
 397         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 398                 if(!strcmp(calibType, fgkDetectorName[iDet])) {
 399                         aPath.SetPath(Form("%s/*", calibType));
 400                         AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
 401                         break;
 402                 }
 403         }
 404         if(!aPath.IsValid()){
 405                 AliError(Form("Not a valid path or detector: %s", calibType));
 406                 return;
 407         }
 408   }
 409
 410 //  // check that calibType refers to a "valid" detector name
 411 //  Bool_t isDetector = kFALSE;
 412 //  for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 413 //    TString detName = fgkDetectorName[iDet];
 414 //    if(aPath.GetLevel0() == detName) {
 415 //      isDetector = kTRUE;
 416 //      break;
 417 //    }
 418 //  }
 419 //
 420 //  if(!isDetector) {
 421 //      AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
 422 //      return;
 423 //  }
 424
 425   TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
 426   if (obj) fSpecCDBUri.Remove(obj);
 427   fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
 428
 429 }
 430
 431
 432
 433
 434 //_____________________________________________________________________________
 435 Bool_t AliReconstruction::SetRunNumber()
 436 {
 437   // The method is called in Run() in order
 438   // to set a correct run number.
 439   // In case of raw data reconstruction the
 440   // run number is taken from the raw data header
 441
 442   if(AliCDBManager::Instance()->GetRun() < 0) {
 443     if (!fRunLoader) {
 444       AliError("No run loader is found !");
 445       return kFALSE;
 446     }
 447     // read run number from gAlice
 448     if(fRunLoader->GetAliRun())
 449       AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
 450     else {
 451       if(fRawReader) {
 452         if(fRawReader->NextEvent()) {
 453           AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
 454           fRawReader->RewindEvents();
 455         }
 456         else {
 457           AliError("No raw-data events found !");
 458           return kFALSE;
 459         }
 460       }
 461       else {
 462         AliError("Neither gAlice nor RawReader objects are found !");
 463         return kFALSE;
 464       }
 465     }
 466     AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
 467   }
 468   return kTRUE;
 469 }
 470
 471 //_____________________________________________________________________________
 472 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
 473 {
 474   // Read the alignment objects from CDB.
 475   // Each detector is supposed to have the
 476   // alignment objects in DET/Align/Data CDB path.
 477   // All the detector objects are then collected,
 478   // sorted by geometry level (starting from ALIC) and
 479   // then applied to the TGeo geometry.
 480   // Finally an overlaps check is performed.
 481
 482   // Load alignment data from CDB and fill fAlignObjArray
 483   if(fLoadAlignFromCDB){
 484
 485     TString detStr = detectors;
 486     TString loadAlObjsListOfDets = "";
 487
 488     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 489       if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 490       loadAlObjsListOfDets += fgkDetectorName[iDet];
 491       loadAlObjsListOfDets += " ";
 492     } // end loop over detectors
 493     loadAlObjsListOfDets.Prepend("GRP "); //add alignment objects for non-sensitive modules
 494     AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
 495   }else{
 496     // Check if the array with alignment objects was
 497     // provided by the user. If yes, apply the objects
 498     // to the present TGeo geometry
 499     if (fAlignObjArray) {
 500       if (gGeoManager && gGeoManager->IsClosed()) {
 501         if (AliGeomManager::ApplyAlignObjsToGeom(*fAlignObjArray) == kFALSE) {
 502           AliError("The misalignment of one or more volumes failed!"
 503                    "Compare the list of simulated detectors and the list of detector alignment data!");
 504           return kFALSE;
 505         }
 506       }
 507       else {
 508         AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
 509         return kFALSE;
 510       }
 511     }
 512   }
 513
 514   delete fAlignObjArray; fAlignObjArray=0;
 515
 516   return kTRUE;
 517 }
 518
 519 //_____________________________________________________________________________
 520 void AliReconstruction::SetGAliceFile(const char* fileName)
 521 {
 522 // set the name of the galice file
 523
 524   fGAliceFileName = fileName;
 525 }
 526
 527 //_____________________________________________________________________________
 528 void AliReconstruction::SetOption(const char* detector, const char* option)
 529 {
 530 // set options for the reconstruction of a detector
 531
 532   TObject* obj = fOptions.FindObject(detector);
 533   if (obj) fOptions.Remove(obj);
 534   fOptions.Add(new TNamed(detector, option));
 535 }
 536
 537
 538 //_____________________________________________________________________________
 539 Bool_t AliReconstruction::Run(const char* input)
 540 {
 541 // run the reconstruction
 542
 543   AliCodeTimerAuto("")
 544
 545   // set the input
 546   if (!input) input = fInput.Data();
 547   TString fileName(input);
 548   if (fileName.EndsWith("/")) {
 549     fRawReader = new AliRawReaderFile(fileName);
 550   } else if (fileName.EndsWith(".root")) {
 551     fRawReader = new AliRawReaderRoot(fileName);
 552   } else if (!fileName.IsNull()) {
 553     fRawReader = new AliRawReaderDate(fileName);
 554     fRawReader->SelectEvents(7);
 555   }
 556   if (!fEquipIdMap.IsNull() && fRawReader)
 557     fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
 558
 559
 560   // get the run loader
 561   if (!InitRunLoader()) return kFALSE;
 562
 563   // Initialize the CDB storage
 564   InitCDBStorage();
 565
 566   // Set run number in CDBManager (if it is not already set by the user)
 567   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
 568
 569   // Import ideal TGeo geometry and apply misalignment
 570   if (!gGeoManager) {
 571     TString geom(gSystem->DirName(fGAliceFileName));
 572     geom += "/geometry.root";
 573     AliGeomManager::LoadGeometry(geom.Data());
 574     if (!gGeoManager) if (fStopOnError) return kFALSE;
 575   }
 576
 577   if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
 578
 579   // local reconstruction
 580   if (!fRunLocalReconstruction.IsNull()) {
 581     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
 582       if (fStopOnError) {CleanUp(); return kFALSE;}
 583     }
 584   }
 585 //  if (!fRunVertexFinder && fRunTracking.IsNull() &&
 586 //      fFillESD.IsNull()) return kTRUE;
 587
 588   // get vertexer
 589   if (fRunVertexFinder && !CreateVertexer()) {
 590     if (fStopOnError) {
 591       CleanUp();
 592       return kFALSE;
 593     }
 594   }
 595
 596   // get trackers
 597   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
 598     if (fStopOnError) {
 599       CleanUp();
 600       return kFALSE;
 601     }
 602   }
 603
 604   // get the possibly already existing ESD file and tree
 605   AliESDEvent* esd = new AliESDEvent(); AliESDEvent* hltesd = new AliESDEvent();
 606   TFile* fileOld = NULL;
 607   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
 608   if (!gSystem->AccessPathName("AliESDs.root")){
 609     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
 610     fileOld = TFile::Open("AliESDs.old.root");
 611     if (fileOld && fileOld->IsOpen()) {
 612       treeOld = (TTree*) fileOld->Get("esdTree");
 613       if (treeOld)esd->ReadFromTree(treeOld);
 614       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
 615       if (hlttreeOld)   hltesd->ReadFromTree(hlttreeOld);
 616     }
 617   }
 618
 619   // create the ESD output file and tree
 620   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 621   file->SetCompressionLevel(2);
 622   if (!file->IsOpen()) {
 623     AliError("opening AliESDs.root failed");
 624     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 625   }
 626
 627   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 628   esd = new AliESDEvent();
 629   esd->CreateStdContent();
 630   esd->WriteToTree(tree);
 631
 632   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
 633   hltesd = new AliESDEvent();
 634   hltesd->CreateStdContent();
 635   hltesd->WriteToTree(hlttree);
 636
 637   /* CKB Why?
 638   delete esd; delete hltesd;
 639   esd = NULL; hltesd = NULL;
 640   */
 641   // create the branch with ESD additions
 642
 643
 644
 645   AliESDfriend *esdf = 0;
 646   if (fWriteESDfriend) {
 647     esdf = new AliESDfriend();
 648     TBranch *br=tree->Branch("ESDfriend.","AliESDfriend", &esdf);
 649     br->SetFile("AliESDfriends.root");
 650     esd->AddObject(esdf);
 651   }
 652
 653
 654   // Get the diamond profile from OCDB
 655   AliCDBEntry* entry = AliCDBManager::Instance()
 656         ->Get("GRP/Calib/MeanVertex");
 657
 658   if(entry) {
 659         fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
 660   } else {
 661         AliError("No diamond profile found in OCDB!");
 662   }
 663
 664   AliVertexerTracks tVertexer(AliTracker::GetBz());
 665   if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
 666
 667   // loop over events
 668
 669   if (fRawReader) fRawReader->RewindEvents();
 670    TString detStr(fFillESD) ;
 671
 672   ProcInfo_t ProcInfo;
 673   gSystem->GetProcInfo(&ProcInfo);
 674   AliInfo(Form("Current memory usage %d %d", ProcInfo.fMemResident, ProcInfo.fMemVirtual));
 675   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 676     if (fRawReader) fRawReader->NextEvent();
 677     if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
 678       // copy old ESD to the new one
 679       if (treeOld) {
 680         esd->ReadFromTree(treeOld);
 681         treeOld->GetEntry(iEvent);
 682       }
 683       tree->Fill();
 684       if (hlttreeOld) {
 685         esd->ReadFromTree(hlttreeOld);
 686         hlttreeOld->GetEntry(iEvent);
 687       }
 688       hlttree->Fill();
 689       continue;
 690     }
 691
 692     AliInfo(Form("processing event %d", iEvent));
 693     fRunLoader->GetEvent(iEvent);
 694
 695     char aFileName[256];
 696     sprintf(aFileName, "ESD_%d.%d_final.root",
 697             fRunLoader->GetHeader()->GetRun(),
 698             fRunLoader->GetHeader()->GetEventNrInRun());
 699     if (!gSystem->AccessPathName(aFileName)) continue;
 700
 701     // local reconstruction
 702     if (!fRunLocalReconstruction.IsNull()) {
 703       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
 704         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 705       }
 706     }
 707
 708
 709     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 710     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 711     esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
 712     hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
 713
 714     // Set magnetic field from the tracker
 715     esd->SetMagneticField(AliTracker::GetBz());
 716     hltesd->SetMagneticField(AliTracker::GetBz());
 717
 718
 719
 720     // Fill raw-data error log into the ESD
 721     if (fRawReader) FillRawDataErrorLog(iEvent,esd);
 722
 723     // vertex finder
 724     if (fRunVertexFinder) {
 725       if (!ReadESD(esd, "vertex")) {
 726         if (!RunVertexFinder(esd)) {
 727           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 728         }
 729         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 730       }
 731     }
 732
 733     // HLT tracking
 734     if (!fRunTracking.IsNull()) {
 735       if (fRunHLTTracking) {
 736         hltesd->SetVertex(esd->GetVertex());
 737         if (!RunHLTTracking(hltesd)) {
 738           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 739         }
 740       }
 741     }
 742
 743     // Muon tracking
 744     if (!fRunTracking.IsNull()) {
 745       if (fRunMuonTracking) {
 746         if (!RunMuonTracking(esd)) {
 747           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 748         }
 749       }
 750     }
 751
 752     // barrel tracking
 753     if (!fRunTracking.IsNull()) {
 754       if (!ReadESD(esd, "tracking")) {
 755         if (!RunTracking(esd)) {
 756           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 757         }
 758         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 759       }
 760     }
 761
 762     // fill ESD
 763     if (!fFillESD.IsNull()) {
 764       if (!FillESD(esd, fFillESD)) {
 765         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 766       }
 767     }
 768
 769     if (!fFillESD.IsNull())
 770     RunQualAss(fFillESD.Data(), esd);
 771
 772     // fill Event header information from the RawEventHeader
 773     if (fRawReader){FillRawEventHeaderESD(esd);}
 774
 775     // combined PID
 776     AliESDpid::MakePID(esd);
 777     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 778
 779     if (fFillTriggerESD) {
 780       if (!ReadESD(esd, "trigger")) {
 781         if (!FillTriggerESD(esd)) {
 782           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 783         }
 784         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
 785       }
 786     }
 787
 788     //Try to improve the reconstructed primary vertex position using the tracks
 789     AliESDVertex *pvtx=0;
 790     Bool_t dovertex=kTRUE;
 791     TObject* obj = fOptions.FindObject("ITS");
 792     if (obj) {
 793       TString optITS = obj->GetTitle();
 794       if (optITS.Contains("cosmics") || optITS.Contains("COSMICS"))
 795         dovertex=kFALSE;
 796     }
 797     if(dovertex) pvtx=tVertexer.FindPrimaryVertex(esd);
 798     if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
 799
 800     if (pvtx)
 801     if (pvtx->GetStatus()) {
 802        // Store the improved primary vertex
 803        esd->SetPrimaryVertex(pvtx);
 804        // Propagate the tracks to the DCA to the improved primary vertex
 805        Double_t somethingbig = 777.;
 806        Double_t bz = esd->GetMagneticField();
 807        Int_t nt=esd->GetNumberOfTracks();
 808        while (nt--) {
 809          AliESDtrack *t = esd->GetTrack(nt);
 810          t->RelateToVertex(pvtx, bz, somethingbig);
 811        }
 812     }
 813
 814     if (fRunV0Finder) {
 815        // V0 finding
 816        AliV0vertexer vtxer;
 817        vtxer.Tracks2V0vertices(esd);
 818
 819        if (fRunCascadeFinder) {
 820           // Cascade finding
 821           AliCascadeVertexer cvtxer;
 822           cvtxer.V0sTracks2CascadeVertices(esd);
 823        }
 824     }
 825
 826     // write ESD
 827     if (fCleanESD) CleanESD(esd);
 828     if (fWriteESDfriend) {
 829       esdf->~AliESDfriend();
 830       new (esdf) AliESDfriend(); // Reset...
 831       esd->GetESDfriend(esdf);
 832     }
 833     tree->Fill();
 834
 835     // write HLT ESD
 836     hlttree->Fill();
 837
 838     if (fCheckPointLevel > 0)  WriteESD(esd, "final");
 839     esd->Reset();
 840     hltesd->Reset();
 841     if (fWriteESDfriend) {
 842       esdf->~AliESDfriend();
 843       new (esdf) AliESDfriend(); // Reset...
 844     }
 845     // esdf->Reset();
 846     // delete esdf; esdf = 0;
 847   // ESD QA
 848
 849     gSystem->GetProcInfo(&ProcInfo);
 850     AliInfo(Form("Event %d -> Current memory usage %d %d",iEvent, ProcInfo.fMemResident, ProcInfo.fMemVirtual));
 851   }
 852
 853   detStr = fFillESD ;
 854   // write quality assurance ESDs data (one entry for all events)
 855   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 856         if (!IsSelected(fgkDetectorName[iDet], detStr))
 857                 continue;
 858     AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
 859     if (!qadm) continue;
 860     qadm->Finish(AliQualAss::kRECPOINTS);
 861     qadm->Finish(AliQualAss::kESDS) ;
 862   }
 863
 864   tree->GetUserInfo()->Add(esd);
 865   hlttree->GetUserInfo()->Add(hltesd);
 866
 867
 868
 869   if(fESDPar.Contains("ESD.par")){
 870     AliInfo("Attaching ESD.par to Tree");
 871     TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
 872     tree->GetUserInfo()->Add(fn);
 873   }
 874
 875
 876   file->cd();
 877   if (fWriteESDfriend)
 878     tree->SetBranchStatus("ESDfriend*",0);
 879   tree->Write();
 880   hlttree->Write();
 881
 882   if (fWriteAOD) {
 883     TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
 884     ESDFile2AODFile(file, aodFile);
 885     aodFile->Close();
 886   }
 887
 888   gROOT->cd();
 889   CleanUp(file, fileOld);
 890
 891   // Create tags for the events in the ESD tree (the ESD tree is always present)
 892   // In case of empty events the tags will contain dummy values
 893   AliESDTagCreator *esdtagCreator = new AliESDTagCreator();
 894   esdtagCreator->CreateESDTags(fFirstEvent,fLastEvent);
 895   if (fWriteAOD) {
 896     AliAODTagCreator *aodtagCreator = new AliAODTagCreator();
 897     aodtagCreator->CreateAODTags(fFirstEvent,fLastEvent);
 898   }
 899
 900   return kTRUE;
 901 }
 902
 903
 904 //_____________________________________________________________________________
 905 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 906 {
 907 // run the local reconstruction
 908
 909   AliCodeTimerAuto("")
 910
 911  //  AliCDBManager* man = AliCDBManager::Instance();
 912 //   Bool_t origCache = man->GetCacheFlag();
 913
 914 //   TString detStr = detectors;
 915 //   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 916 //     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 917 //     AliReconstructor* reconstructor = GetReconstructor(iDet);
 918 //     if (!reconstructor) continue;
 919 //     if (reconstructor->HasLocalReconstruction()) continue;
 920
 921 //     AliCodeTimerStart(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 922 //     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 923
 924 //     AliCodeTimerStart(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
 925 //     AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
 926
 927 //     man->SetCacheFlag(kTRUE);
 928 //     TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
 929 //     man->GetAll(calibPath); // entries are cached!
 930
 931 //     AliCodeTimerStop(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
 932
 933 //     if (fRawReader) {
 934 //       fRawReader->RewindEvents();
 935 //       reconstructor->Reconstruct(fRunLoader, fRawReader);
 936 //     } else {
 937 //       reconstructor->Reconstruct(fRunLoader);
 938 //     }
 939
 940 //      AliCodeTimerStop(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 941
 942 //     // unload calibration data
 943 //     man->UnloadFromCache(calibPath);
 944 //     //man->ClearCache();
 945 //   }
 946
 947 //   man->SetCacheFlag(origCache);
 948
 949 //   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 950 //     AliError(Form("the following detectors were not found: %s",
 951 //                   detStr.Data()));
 952 //     if (fStopOnError) return kFALSE;
 953 //   }
 954
 955   return kTRUE;
 956 }
 957
 958 //_____________________________________________________________________________
 959 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 960 {
 961 // run the local reconstruction
 962
 963   AliCodeTimerAuto("")
 964
 965   TString detStr = detectors;
 966   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 967     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 968     AliReconstructor* reconstructor = GetReconstructor(iDet);
 969     if (!reconstructor) continue;
 970     AliLoader* loader = fLoader[iDet];
 971     if (!loader) {
 972       AliWarning(Form("No loader is defined for %s!",fgkDetectorName[iDet]));
 973       continue;
 974     }
 975
 976     // conversion of digits
 977     if (fRawReader && reconstructor->HasDigitConversion()) {
 978       AliInfo(Form("converting raw data digits into root objects for %s",
 979                    fgkDetectorName[iDet]));
 980       AliCodeTimerAuto(Form("converting raw data digits into root objects for %s",
 981                             fgkDetectorName[iDet]));
 982       loader->LoadDigits("update");
 983       loader->CleanDigits();
 984       loader->MakeDigitsContainer();
 985       TTree* digitsTree = loader->TreeD();
 986       reconstructor->ConvertDigits(fRawReader, digitsTree);
 987       loader->WriteDigits("OVERWRITE");
 988       loader->UnloadDigits();
 989     }
 990
 991     // local reconstruction
 992     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 993     AliCodeTimerAuto(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 994     loader->LoadRecPoints("update");
 995     loader->CleanRecPoints();
 996     loader->MakeRecPointsContainer();
 997     TTree* clustersTree = loader->TreeR();
 998     if (fRawReader && !reconstructor->HasDigitConversion()) {
 999       reconstructor->Reconstruct(fRawReader, clustersTree);
1000     } else {
1001       loader->LoadDigits("read");
1002       TTree* digitsTree = loader->TreeD();
1003       if (!digitsTree) {
1004         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
1005         if (fStopOnError) return kFALSE;
1006       } else {
1007         reconstructor->Reconstruct(digitsTree, clustersTree);
1008       }
1009       loader->UnloadDigits();
1010     }
1011
1012     AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
1013     if (qadm) {
1014       AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1015       AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1016       qadm->Exec(AliQualAss::kRECPOINTS, clustersTree) ;
1017       AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1018     }
1019
1020     loader->WriteRecPoints("OVERWRITE");
1021     loader->UnloadRecPoints();
1022   }
1023
1024   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1025     AliError(Form("the following detectors were not found: %s",
1026                   detStr.Data()));
1027     if (fStopOnError) return kFALSE;
1028   }
1029
1030   return kTRUE;
1031 }
1032
1033 //_____________________________________________________________________________
1034 Bool_t AliReconstruction::RunVertexFinder(AliESDEvent*& esd)
1035 {
1036 // run the barrel tracking
1037
1038   AliCodeTimerAuto("")
1039
1040   AliESDVertex* vertex = NULL;
1041   Double_t vtxPos[3] = {0, 0, 0};
1042   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
1043   TArrayF mcVertex(3);
1044   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
1045     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
1046     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
1047   }
1048
1049   if (fVertexer) {
1050     if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
1051     AliInfo("running the ITS vertex finder");
1052     if (fLoader[0]) fLoader[0]->LoadRecPoints();
1053     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
1054     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
1055     if(!vertex){
1056       AliWarning("Vertex not found");
1057       vertex = new AliESDVertex();
1058       vertex->SetName("default");
1059     }
1060     else {
1061       vertex->SetName("reconstructed");
1062     }
1063
1064   } else {
1065     AliInfo("getting the primary vertex from MC");
1066     vertex = new AliESDVertex(vtxPos, vtxErr);
1067   }
1068
1069   if (vertex) {
1070     vertex->GetXYZ(vtxPos);
1071     vertex->GetSigmaXYZ(vtxErr);
1072   } else {
1073     AliWarning("no vertex reconstructed");
1074     vertex = new AliESDVertex(vtxPos, vtxErr);
1075   }
1076   esd->SetVertex(vertex);
1077   // if SPD multiplicity has been determined, it is stored in the ESD
1078   AliMultiplicity *mult = fVertexer->GetMultiplicity();
1079   if(mult)esd->SetMultiplicity(mult);
1080
1081   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1082     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1083   }
1084   delete vertex;
1085
1086   return kTRUE;
1087 }
1088
1089 //_____________________________________________________________________________
1090 Bool_t AliReconstruction::RunHLTTracking(AliESDEvent*& esd)
1091 {
1092 // run the HLT barrel tracking
1093
1094   AliCodeTimerAuto("")
1095
1096   if (!fRunLoader) {
1097     AliError("Missing runLoader!");
1098     return kFALSE;
1099   }
1100
1101   AliInfo("running HLT tracking");
1102
1103   // Get a pointer to the HLT reconstructor
1104   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1105   if (!reconstructor) return kFALSE;
1106
1107   // TPC + ITS
1108   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1109     TString detName = fgkDetectorName[iDet];
1110     AliDebug(1, Form("%s HLT tracking", detName.Data()));
1111     reconstructor->SetOption(detName.Data());
1112     AliTracker *tracker = reconstructor->CreateTracker();
1113     if (!tracker) {
1114       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1115       if (fStopOnError) return kFALSE;
1116       continue;
1117     }
1118     Double_t vtxPos[3];
1119     Double_t vtxErr[3]={0.005,0.005,0.010};
1120     const AliESDVertex *vertex = esd->GetVertex();
1121     vertex->GetXYZ(vtxPos);
1122     tracker->SetVertex(vtxPos,vtxErr);
1123     if(iDet != 1) {
1124       fLoader[iDet]->LoadRecPoints("read");
1125       TTree* tree = fLoader[iDet]->TreeR();
1126       if (!tree) {
1127         AliError(Form("Can't get the %s cluster tree", detName.Data()));
1128         return kFALSE;
1129       }
1130       tracker->LoadClusters(tree);
1131     }
1132     if (tracker->Clusters2Tracks(esd) != 0) {
1133       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1134       return kFALSE;
1135     }
1136     if(iDet != 1) {
1137       tracker->UnloadClusters();
1138     }
1139     delete tracker;
1140   }
1141
1142   return kTRUE;
1143 }
1144
1145 //_____________________________________________________________________________
1146 Bool_t AliReconstruction::RunMuonTracking(AliESDEvent*& esd)
1147 {
1148 // run the muon spectrometer tracking
1149
1150   AliCodeTimerAuto("")
1151
1152   if (!fRunLoader) {
1153     AliError("Missing runLoader!");
1154     return kFALSE;
1155   }
1156   Int_t iDet = 7; // for MUON
1157
1158   AliInfo("is running...");
1159
1160   // Get a pointer to the MUON reconstructor
1161   AliReconstructor *reconstructor = GetReconstructor(iDet);
1162   if (!reconstructor) return kFALSE;
1163
1164
1165   TString detName = fgkDetectorName[iDet];
1166   AliDebug(1, Form("%s tracking", detName.Data()));
1167   AliTracker *tracker =  reconstructor->CreateTracker();
1168   if (!tracker) {
1169     AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1170     return kFALSE;
1171   }
1172
1173   // create Tracks
1174   fLoader[iDet]->LoadTracks("update");
1175   fLoader[iDet]->CleanTracks();
1176   fLoader[iDet]->MakeTracksContainer();
1177
1178   // read RecPoints
1179   fLoader[iDet]->LoadRecPoints("read");
1180   tracker->LoadClusters(fLoader[iDet]->TreeR());
1181
1182   Int_t rv = tracker->Clusters2Tracks(esd);
1183
1184   fLoader[iDet]->UnloadRecPoints();
1185
1186   if ( rv )
1187   {
1188     AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1189     return kFALSE;
1190   }
1191
1192   tracker->UnloadClusters();
1193
1194   fLoader[iDet]->UnloadRecPoints();
1195
1196   fLoader[iDet]->WriteTracks("OVERWRITE");
1197   fLoader[iDet]->UnloadTracks();
1198
1199   delete tracker;
1200
1201   return kTRUE;
1202 }
1203
1204
1205 //_____________________________________________________________________________
1206 Bool_t AliReconstruction::RunTracking(AliESDEvent*& esd)
1207 {
1208 // run the barrel tracking
1209
1210   AliCodeTimerAuto("")
1211
1212   AliInfo("running tracking");
1213
1214   //Fill the ESD with the T0 info (will be used by the TOF)
1215   if (fReconstructor[11] && fLoader[11]) {
1216     fLoader[11]->LoadRecPoints("READ");
1217     TTree *treeR = fLoader[11]->TreeR();
1218     GetReconstructor(11)->FillESD((TTree *)NULL,treeR,esd);
1219   }
1220
1221   // pass 1: TPC + ITS inwards
1222   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1223     if (!fTracker[iDet]) continue;
1224     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1225
1226     // load clusters
1227     fLoader[iDet]->LoadRecPoints("read");
1228     TTree* tree = fLoader[iDet]->TreeR();
1229     if (!tree) {
1230       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1231       return kFALSE;
1232     }
1233     fTracker[iDet]->LoadClusters(tree);
1234
1235     // run tracking
1236     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1237       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1238       return kFALSE;
1239     }
1240     if (fCheckPointLevel > 1) {
1241       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1242     }
1243     // preliminary PID in TPC needed by the ITS tracker
1244     if (iDet == 1) {
1245       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1246       AliESDpid::MakePID(esd);
1247     }
1248   }
1249
1250   // pass 2: ALL backwards
1251   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1252     if (!fTracker[iDet]) continue;
1253     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1254
1255     // load clusters
1256     if (iDet > 1) {     // all except ITS, TPC
1257       TTree* tree = NULL;
1258       fLoader[iDet]->LoadRecPoints("read");
1259       tree = fLoader[iDet]->TreeR();
1260       if (!tree) {
1261         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1262         return kFALSE;
1263       }
1264       fTracker[iDet]->LoadClusters(tree);
1265     }
1266
1267     // run tracking
1268     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1269       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1270       //      return kFALSE;
1271     }
1272     if (fCheckPointLevel > 1) {
1273       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1274     }
1275
1276     // unload clusters
1277     if (iDet > 2) {     // all except ITS, TPC, TRD
1278       fTracker[iDet]->UnloadClusters();
1279       fLoader[iDet]->UnloadRecPoints();
1280     }
1281     // updated PID in TPC needed by the ITS tracker -MI
1282     if (iDet == 1) {
1283       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1284       AliESDpid::MakePID(esd);
1285     }
1286   }
1287
1288   // write space-points to the ESD in case alignment data output
1289   // is switched on
1290   if (fWriteAlignmentData)
1291     WriteAlignmentData(esd);
1292
1293   // pass 3: TRD + TPC + ITS refit inwards
1294   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1295     if (!fTracker[iDet]) continue;
1296     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1297
1298     // run tracking
1299     if (fTracker[iDet]->RefitInward(esd) != 0) {
1300       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1301       //      return kFALSE;
1302     }
1303     if (fCheckPointLevel > 1) {
1304       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1305     }
1306
1307     // unload clusters
1308     fTracker[iDet]->UnloadClusters();
1309     fLoader[iDet]->UnloadRecPoints();
1310   }
1311   //
1312   // Propagate track to the vertex - if not done by ITS
1313   //
1314   Int_t ntracks = esd->GetNumberOfTracks();
1315   for (Int_t itrack=0; itrack<ntracks; itrack++){
1316     const Double_t kRadius  = 3;   // beam pipe radius
1317     const Double_t kMaxStep = 5;   // max step
1318     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1319     Double_t       fieldZ   = AliTracker::GetBz();  //
1320     AliESDtrack * track = esd->GetTrack(itrack);
1321     if (!track) continue;
1322     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1323    AliTracker::PropagateTrackTo(track,kRadius,track->GetMass(),kMaxStep,kTRUE);
1324     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1325   }
1326
1327   return kTRUE;
1328 }
1329
1330 //_____________________________________________________________________________
1331 Bool_t AliReconstruction::CleanESD(AliESDEvent *esd){
1332   //
1333   // Remove the data which are not needed for the physics analysis.
1334   //
1335
1336   AliInfo("Cleaning the ESD...");
1337
1338   const AliESDVertex *vertex=esd->GetVertex();
1339   Double_t vz=vertex->GetZv();
1340
1341   Int_t nTracks=esd->GetNumberOfTracks();
1342   for (Int_t i=0; i<nTracks; i++) {
1343     AliESDtrack *track=esd->GetTrack(i);
1344
1345     Float_t xy,z; track->GetImpactParameters(xy,z);
1346     if (TMath::Abs(xy) < 50.)    continue;
1347     if (vertex->GetStatus())
1348       if (TMath::Abs(vz-z) < 5.) continue;
1349
1350     esd->RemoveTrack(i);
1351   }
1352
1353   return kTRUE;
1354 }
1355
1356 //_____________________________________________________________________________
1357 Bool_t AliReconstruction::FillESD(AliESDEvent*& esd, const TString& detectors)
1358 {
1359 // fill the event summary data
1360
1361   AliCodeTimerAuto("")
1362
1363   TString detStr = detectors;
1364   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1365     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1366     AliReconstructor* reconstructor = GetReconstructor(iDet);
1367     if (!reconstructor) continue;
1368
1369     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1370       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1371       TTree* clustersTree = NULL;
1372       if (fLoader[iDet]) {
1373         fLoader[iDet]->LoadRecPoints("read");
1374         clustersTree = fLoader[iDet]->TreeR();
1375         if (!clustersTree) {
1376           AliError(Form("Can't get the %s clusters tree",
1377                         fgkDetectorName[iDet]));
1378           if (fStopOnError) return kFALSE;
1379         }
1380       }
1381       if (fRawReader && !reconstructor->HasDigitConversion()) {
1382         reconstructor->FillESD(fRawReader, clustersTree, esd);
1383       } else {
1384         TTree* digitsTree = NULL;
1385         if (fLoader[iDet]) {
1386           fLoader[iDet]->LoadDigits("read");
1387           digitsTree = fLoader[iDet]->TreeD();
1388           if (!digitsTree) {
1389             AliError(Form("Can't get the %s digits tree",
1390                           fgkDetectorName[iDet]));
1391             if (fStopOnError) return kFALSE;
1392           }
1393         }
1394         reconstructor->FillESD(digitsTree, clustersTree, esd);
1395         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1396       }
1397       if (fLoader[iDet]) {
1398         fLoader[iDet]->UnloadRecPoints();
1399       }
1400
1401       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1402     }
1403   }
1404
1405   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1406     AliError(Form("the following detectors were not found: %s",
1407                   detStr.Data()));
1408     if (fStopOnError) return kFALSE;
1409   }
1410
1411   return kTRUE;
1412 }
1413
1414 //_____________________________________________________________________________
1415 Bool_t AliReconstruction::FillTriggerESD(AliESDEvent*& esd)
1416 {
1417   // Reads the trigger decision which is
1418   // stored in Trigger.root file and fills
1419   // the corresponding esd entries
1420
1421   AliCodeTimerAuto("")
1422
1423   AliInfo("Filling trigger information into the ESD");
1424
1425   if (fRawReader) {
1426     AliCTPRawStream input(fRawReader);
1427     if (!input.Next()) {
1428       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1429       return kFALSE;
1430     }
1431     esd->SetTriggerMask(input.GetClassMask());
1432     esd->SetTriggerCluster(input.GetClusterMask());
1433   }
1434   else {
1435     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1436     if (runloader) {
1437       if (!runloader->LoadTrigger()) {
1438         AliCentralTrigger *aCTP = runloader->GetTrigger();
1439         esd->SetTriggerMask(aCTP->GetClassMask());
1440         esd->SetTriggerCluster(aCTP->GetClusterMask());
1441       }
1442       else {
1443         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1444         return kFALSE;
1445       }
1446     }
1447     else {
1448       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1449       return kFALSE;
1450     }
1451   }
1452
1453   return kTRUE;
1454 }
1455
1456
1457
1458
1459
1460 //_____________________________________________________________________________
1461 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESDEvent*& esd)
1462 {
1463   //
1464   // Filling information from RawReader Header
1465   //
1466
1467   AliInfo("Filling information from RawReader Header");
1468   esd->SetBunchCrossNumber(0);
1469   esd->SetOrbitNumber(0);
1470   esd->SetPeriodNumber(0);
1471   esd->SetTimeStamp(0);
1472   esd->SetEventType(0);
1473   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1474   if (eventHeader){
1475
1476     const UInt_t *id = eventHeader->GetP("Id");
1477     esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1478     esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1479     esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1480
1481     esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1482     esd->SetEventType((eventHeader->Get("Type")));
1483   }
1484
1485   return kTRUE;
1486 }
1487
1488
1489 //_____________________________________________________________________________
1490 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1491 {
1492 // check whether detName is contained in detectors
1493 // if yes, it is removed from detectors
1494
1495   // check if all detectors are selected
1496   if ((detectors.CompareTo("ALL") == 0) ||
1497       detectors.BeginsWith("ALL ") ||
1498       detectors.EndsWith(" ALL") ||
1499       detectors.Contains(" ALL ")) {
1500     detectors = "ALL";
1501     return kTRUE;
1502   }
1503
1504   // search for the given detector
1505   Bool_t result = kFALSE;
1506   if ((detectors.CompareTo(detName) == 0) ||
1507       detectors.BeginsWith(detName+" ") ||
1508       detectors.EndsWith(" "+detName) ||
1509       detectors.Contains(" "+detName+" ")) {
1510     detectors.ReplaceAll(detName, "");
1511     result = kTRUE;
1512   }
1513
1514   // clean up the detectors string
1515   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1516   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1517   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1518
1519   return result;
1520 }
1521
1522 //_____________________________________________________________________________
1523 Bool_t AliReconstruction::InitRunLoader()
1524 {
1525 // get or create the run loader
1526
1527   if (gAlice) delete gAlice;
1528   gAlice = NULL;
1529
1530   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1531     // load all base libraries to get the loader classes
1532     TString libs = gSystem->GetLibraries();
1533     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1534       TString detName = fgkDetectorName[iDet];
1535       if (detName == "HLT") continue;
1536       if (libs.Contains("lib" + detName + "base.so")) continue;
1537       gSystem->Load("lib" + detName + "base.so");
1538     }
1539     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1540     if (!fRunLoader) {
1541       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1542       CleanUp();
1543       return kFALSE;
1544     }
1545     fRunLoader->CdGAFile();
1546     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1547       if (fRunLoader->LoadgAlice() == 0) {
1548         gAlice = fRunLoader->GetAliRun();
1549         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1550       }
1551     }
1552     if (!gAlice && !fRawReader) {
1553       AliError(Form("no gAlice object found in file %s",
1554                     fGAliceFileName.Data()));
1555       CleanUp();
1556       return kFALSE;
1557     }
1558
1559     //PH This is a temporary fix to give access to the kinematics
1560     //PH that is needed for the labels of ITS clusters
1561     fRunLoader->LoadKinematics();
1562
1563   } else {               // galice.root does not exist
1564     if (!fRawReader) {
1565       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1566       CleanUp();
1567       return kFALSE;
1568     }
1569     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1570                                     AliConfig::GetDefaultEventFolderName(),
1571                                     "recreate");
1572     if (!fRunLoader) {
1573       AliError(Form("could not create run loader in file %s",
1574                     fGAliceFileName.Data()));
1575       CleanUp();
1576       return kFALSE;
1577     }
1578     fRunLoader->MakeTree("E");
1579     Int_t iEvent = 0;
1580     while (fRawReader->NextEvent()) {
1581       fRunLoader->SetEventNumber(iEvent);
1582       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1583                                      iEvent, iEvent);
1584       fRunLoader->MakeTree("H");
1585       fRunLoader->TreeE()->Fill();
1586       iEvent++;
1587     }
1588     fRawReader->RewindEvents();
1589     if (fNumberOfEventsPerFile > 0)
1590       fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1591     else
1592       fRunLoader->SetNumberOfEventsPerFile(iEvent);
1593     fRunLoader->WriteHeader("OVERWRITE");
1594     fRunLoader->CdGAFile();
1595     fRunLoader->Write(0, TObject::kOverwrite);
1596 //    AliTracker::SetFieldMap(???);
1597   }
1598
1599   return kTRUE;
1600 }
1601
1602 //_____________________________________________________________________________
1603 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1604 {
1605 // get the reconstructor object and the loader for a detector
1606
1607   if (fReconstructor[iDet]) return fReconstructor[iDet];
1608
1609   // load the reconstructor object
1610   TPluginManager* pluginManager = gROOT->GetPluginManager();
1611   TString detName = fgkDetectorName[iDet];
1612   TString recName = "Ali" + detName + "Reconstructor";
1613   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1614
1615   if (detName == "HLT") {
1616     if (!gROOT->GetClass("AliLevel3")) {
1617       gSystem->Load("libAliHLTSrc.so");
1618       gSystem->Load("libAliHLTMisc.so");
1619       gSystem->Load("libAliHLTHough.so");
1620       gSystem->Load("libAliHLTComp.so");
1621     }
1622   }
1623
1624   AliReconstructor* reconstructor = NULL;
1625   // first check if a plugin is defined for the reconstructor
1626   TPluginHandler* pluginHandler =
1627     pluginManager->FindHandler("AliReconstructor", detName);
1628   // if not, add a plugin for it
1629   if (!pluginHandler) {
1630     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1631     TString libs = gSystem->GetLibraries();
1632     if (libs.Contains("lib" + detName + "base.so") ||
1633         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1634       pluginManager->AddHandler("AliReconstructor", detName,
1635                                 recName, detName + "rec", recName + "()");
1636     } else {
1637       pluginManager->AddHandler("AliReconstructor", detName,
1638                                 recName, detName, recName + "()");
1639     }
1640     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1641   }
1642   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1643     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1644   }
1645   if (reconstructor) {
1646     TObject* obj = fOptions.FindObject(detName.Data());
1647     if (obj) reconstructor->SetOption(obj->GetTitle());
1648     reconstructor->Init();
1649     fReconstructor[iDet] = reconstructor;
1650   }
1651
1652   // get or create the loader
1653   if (detName != "HLT") {
1654     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1655     if (!fLoader[iDet]) {
1656       AliConfig::Instance()
1657         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1658                                 detName, detName);
1659       // first check if a plugin is defined for the loader
1660       pluginHandler =
1661         pluginManager->FindHandler("AliLoader", detName);
1662       // if not, add a plugin for it
1663       if (!pluginHandler) {
1664         TString loaderName = "Ali" + detName + "Loader";
1665         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1666         pluginManager->AddHandler("AliLoader", detName,
1667                                   loaderName, detName + "base",
1668                                   loaderName + "(const char*, TFolder*)");
1669         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1670       }
1671       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1672         fLoader[iDet] =
1673           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1674                                                  fRunLoader->GetEventFolder());
1675       }
1676       if (!fLoader[iDet]) {   // use default loader
1677         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1678       }
1679       if (!fLoader[iDet]) {
1680         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1681         if (fStopOnError) return NULL;
1682       } else {
1683         fRunLoader->AddLoader(fLoader[iDet]);
1684         fRunLoader->CdGAFile();
1685         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1686         fRunLoader->Write(0, TObject::kOverwrite);
1687       }
1688     }
1689   }
1690
1691   return reconstructor;
1692 }
1693
1694 //_____________________________________________________________________________
1695 Bool_t AliReconstruction::CreateVertexer()
1696 {
1697 // create the vertexer
1698
1699   fVertexer = NULL;
1700   AliReconstructor* itsReconstructor = GetReconstructor(0);
1701   if (itsReconstructor) {
1702     fVertexer = itsReconstructor->CreateVertexer();
1703   }
1704   if (!fVertexer) {
1705     AliWarning("couldn't create a vertexer for ITS");
1706     if (fStopOnError) return kFALSE;
1707   }
1708
1709   return kTRUE;
1710 }
1711
1712 //_____________________________________________________________________________
1713 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1714 {
1715 // create the trackers
1716
1717   TString detStr = detectors;
1718   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1719     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1720     AliReconstructor* reconstructor = GetReconstructor(iDet);
1721     if (!reconstructor) continue;
1722     TString detName = fgkDetectorName[iDet];
1723     if (detName == "HLT") {
1724       fRunHLTTracking = kTRUE;
1725       continue;
1726     }
1727     if (detName == "MUON") {
1728       fRunMuonTracking = kTRUE;
1729       continue;
1730     }
1731
1732
1733     fTracker[iDet] = reconstructor->CreateTracker();
1734     if (!fTracker[iDet] && (iDet < 7)) {
1735       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1736       if (fStopOnError) return kFALSE;
1737     }
1738   }
1739
1740   return kTRUE;
1741 }
1742
1743 //_____________________________________________________________________________
1744 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1745 {
1746 // delete trackers and the run loader and close and delete the file
1747
1748   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1749     delete fReconstructor[iDet];
1750     fReconstructor[iDet] = NULL;
1751     fLoader[iDet] = NULL;
1752     delete fTracker[iDet];
1753     fTracker[iDet] = NULL;
1754     delete fQualAssDataMaker[iDet];
1755     fQualAssDataMaker[iDet] = NULL;
1756   }
1757   delete fVertexer;
1758   fVertexer = NULL;
1759   delete fDiamondProfile;
1760   fDiamondProfile = NULL;
1761
1762   delete fRunLoader;
1763   fRunLoader = NULL;
1764   delete fRawReader;
1765   fRawReader = NULL;
1766
1767   if (file) {
1768     file->Close();
1769     delete file;
1770   }
1771
1772   if (fileOld) {
1773     fileOld->Close();
1774     delete fileOld;
1775     gSystem->Unlink("AliESDs.old.root");
1776   }
1777 }
1778
1779 //_____________________________________________________________________________
1780
1781 Bool_t AliReconstruction::ReadESD(AliESDEvent*& esd, const char* recStep) const
1782 {
1783 // read the ESD event from a file
1784
1785   if (!esd) return kFALSE;
1786   char fileName[256];
1787   sprintf(fileName, "ESD_%d.%d_%s.root",
1788           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1789   if (gSystem->AccessPathName(fileName)) return kFALSE;
1790
1791   AliInfo(Form("reading ESD from file %s", fileName));
1792   AliDebug(1, Form("reading ESD from file %s", fileName));
1793   TFile* file = TFile::Open(fileName);
1794   if (!file || !file->IsOpen()) {
1795     AliError(Form("opening %s failed", fileName));
1796     delete file;
1797     return kFALSE;
1798   }
1799
1800   gROOT->cd();
1801   delete esd;
1802   esd = (AliESDEvent*) file->Get("ESD");
1803   file->Close();
1804   delete file;
1805   return kTRUE;
1806
1807 }
1808
1809
1810
1811 //_____________________________________________________________________________
1812 void AliReconstruction::WriteESD(AliESDEvent* esd, const char* recStep) const
1813 {
1814 // write the ESD event to a file
1815
1816   if (!esd) return;
1817   char fileName[256];
1818   sprintf(fileName, "ESD_%d.%d_%s.root",
1819           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1820
1821   AliDebug(1, Form("writing ESD to file %s", fileName));
1822   TFile* file = TFile::Open(fileName, "recreate");
1823   if (!file || !file->IsOpen()) {
1824     AliError(Form("opening %s failed", fileName));
1825   } else {
1826     esd->Write("ESD");
1827     file->Close();
1828   }
1829   delete file;
1830 }
1831
1832
1833
1834
1835
1836 //_____________________________________________________________________________
1837 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
1838 {
1839   // write all files from the given esd file to an aod file
1840
1841   // create an AliAOD object
1842   AliAODEvent *aod = new AliAODEvent();
1843   aod->CreateStdContent();
1844
1845   // go to the file
1846   aodFile->cd();
1847
1848   // create the tree
1849   TTree *aodTree = new TTree("aodTree", "AliAOD tree");
1850   aodTree->Branch(aod->GetList());
1851
1852   // connect to ESD
1853   TTree *t = (TTree*) esdFile->Get("esdTree");
1854   AliESDEvent *esd = new AliESDEvent();
1855   esd->ReadFromTree(t);
1856
1857   Int_t nEvents = t->GetEntries();
1858
1859   // set arrays and pointers
1860   Float_t posF[3];
1861   Double_t pos[3];
1862   Double_t p[3];
1863   Double_t covVtx[6];
1864   Double_t covTr[21];
1865   Double_t pid[10];
1866
1867   // loop over events and fill them
1868   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
1869     t->GetEntry(iEvent);
1870
1871     // Multiplicity information needed by the header (to be revised!)
1872     Int_t nTracks   = esd->GetNumberOfTracks();
1873     Int_t nPosTracks = 0;
1874     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
1875       if (esd->GetTrack(iTrack)->Charge()> 0) nPosTracks++;
1876
1877     // Access the header
1878     AliAODHeader *header = aod->GetHeader();
1879
1880     // fill the header
1881     header->SetRunNumber       (esd->GetRunNumber()       );
1882     header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
1883     header->SetOrbitNumber     (esd->GetOrbitNumber()     );
1884     header->SetPeriodNumber    (esd->GetPeriodNumber()    );
1885     header->SetTriggerMask     (esd->GetTriggerMask()     );
1886     header->SetTriggerCluster  (esd->GetTriggerCluster()  );
1887     header->SetEventType       (esd->GetEventType()       );
1888     header->SetMagneticField   (esd->GetMagneticField()   );
1889     header->SetZDCN1Energy     (esd->GetZDCN1Energy()     );
1890     header->SetZDCP1Energy     (esd->GetZDCP1Energy()     );
1891     header->SetZDCN2Energy     (esd->GetZDCN2Energy()     );
1892     header->SetZDCP2Energy     (esd->GetZDCP2Energy()     );
1893     header->SetZDCEMEnergy     (esd->GetZDCEMEnergy()     );
1894     header->SetRefMultiplicity   (nTracks);
1895     header->SetRefMultiplicityPos(nPosTracks);
1896     header->SetRefMultiplicityNeg(nTracks - nPosTracks);
1897     header->SetMuonMagFieldScale(-999.); // FIXME
1898     header->SetCentrality(-999.);        // FIXME
1899
1900     Int_t nV0s      = esd->GetNumberOfV0s();
1901     Int_t nCascades = esd->GetNumberOfCascades();
1902     Int_t nKinks    = esd->GetNumberOfKinks();
1903     Int_t nVertices = nV0s + nCascades + nKinks;
1904
1905     aod->ResetStd(nTracks, nVertices);
1906     AliAODTrack *aodTrack;
1907
1908     // Array to take into account the tracks already added to the AOD
1909     Bool_t * usedTrack = NULL;
1910     if (nTracks>0) {
1911       usedTrack = new Bool_t[nTracks];
1912       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
1913     }
1914     // Array to take into account the V0s already added to the AOD
1915     Bool_t * usedV0 = NULL;
1916     if (nV0s>0) {
1917       usedV0 = new Bool_t[nV0s];
1918       for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
1919     }
1920     // Array to take into account the kinks already added to the AOD
1921     Bool_t * usedKink = NULL;
1922     if (nKinks>0) {
1923       usedKink = new Bool_t[nKinks];
1924       for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
1925     }
1926
1927     // Access to the AOD container of vertices
1928     TClonesArray &vertices = *(aod->GetVertices());
1929     Int_t jVertices=0;
1930
1931     // Access to the AOD container of tracks
1932     TClonesArray &tracks = *(aod->GetTracks());
1933     Int_t jTracks=0;
1934
1935     // Add primary vertex. The primary tracks will be defined
1936     // after the loops on the composite objects (V0, cascades, kinks)
1937     const AliESDVertex *vtx = esd->GetPrimaryVertex();
1938
1939     vtx->GetXYZ(pos); // position
1940     vtx->GetCovMatrix(covVtx); //covariance matrix
1941
1942     AliAODVertex * primary = new(vertices[jVertices++])
1943       AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1944
1945     // Create vertices starting from the most complex objects
1946
1947     // Cascades
1948     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
1949       AliESDcascade *cascade = esd->GetCascade(nCascade);
1950
1951       cascade->GetXYZcascade(pos[0], pos[1], pos[2]);  // Bo: bug correction
1952       cascade->GetPosCovXi(covVtx);
1953
1954       // Add the cascade vertex
1955       AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
1956                                                                         covVtx,
1957                                                                         cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
1958                                                                         primary,
1959                                                                         nCascade,
1960                                                                         AliAODVertex::kCascade);
1961
1962       primary->AddDaughter(vcascade);
1963
1964       // Add the V0 from the cascade. The ESD class have to be optimized...
1965       // Now we have to search for the corresponding V0 in the list of V0s
1966       // using the indeces of the positive and negative tracks
1967
1968       Int_t posFromV0 = cascade->GetPindex();
1969       Int_t negFromV0 = cascade->GetNindex();
1970
1971       AliESDv0 * v0 = 0x0;
1972       Int_t indV0 = -1;
1973
1974       for (Int_t iV0=0; iV0<nV0s; ++iV0) {
1975
1976         v0 = esd->GetV0(iV0);
1977         Int_t posV0 = v0->GetPindex();
1978         Int_t negV0 = v0->GetNindex();
1979
1980         if (posV0==posFromV0 && negV0==negFromV0) {
1981           indV0 = iV0;
1982           break;
1983         }
1984       }
1985
1986       AliAODVertex * vV0FromCascade = 0x0;
1987
1988       if (indV0>-1 && !usedV0[indV0] ) {
1989
1990         // the V0 exists in the array of V0s and is not used
1991
1992         usedV0[indV0] = kTRUE;
1993
1994         v0->GetXYZ(pos[0], pos[1], pos[2]);
1995         v0->GetPosCov(covVtx);
1996
1997         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
1998                                                                  covVtx,
1999                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2000                                                                  vcascade,
2001                                                                  indV0,
2002                                                                  AliAODVertex::kV0);
2003       } else {
2004
2005         // the V0 doesn't exist in the array of V0s or was used
2006         cerr << "Error: event " << iEvent << " cascade " << nCascade
2007              << " The V0 " << indV0
2008              << " doesn't exist in the array of V0s or was used!" << endl;
2009
2010         cascade->GetXYZ(pos[0], pos[1], pos[2]);
2011         cascade->GetPosCov(covVtx);
2012
2013         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2014                                                                  covVtx,
2015                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2016                                                                  vcascade,
2017                                                                  indV0,
2018                                                                  AliAODVertex::kV0);
2019         vcascade->AddDaughter(vV0FromCascade);
2020       }
2021
2022       // Add the positive tracks from the V0
2023
2024       if (! usedTrack[posFromV0]) {
2025
2026         usedTrack[posFromV0] = kTRUE;
2027
2028         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2029         esdTrack->GetPxPyPz(p);
2030         esdTrack->GetXYZ(pos);
2031         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2032         esdTrack->GetESDpid(pid);
2033
2034         vV0FromCascade->AddDaughter(aodTrack =
2035                                     new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2036                                            esdTrack->GetLabel(),
2037                                            p,
2038                                            kTRUE,
2039                                            pos,
2040                                            kFALSE,
2041                                            covTr,
2042                                            (Short_t)esdTrack->Charge(),
2043                                            esdTrack->GetITSClusterMap(),
2044                                            pid,
2045                                            vV0FromCascade,
2046                                            kTRUE,  // check if this is right
2047                                            kFALSE, // check if this is right
2048                                            AliAODTrack::kSecondary)
2049                 );
2050         aodTrack->ConvertAliPIDtoAODPID();
2051       }
2052       else {
2053         cerr << "Error: event " << iEvent << " cascade " << nCascade
2054              << " track " << posFromV0 << " has already been used!" << endl;
2055       }
2056
2057       // Add the negative tracks from the V0
2058
2059       if (!usedTrack[negFromV0]) {
2060
2061         usedTrack[negFromV0] = kTRUE;
2062
2063         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2064         esdTrack->GetPxPyPz(p);
2065         esdTrack->GetXYZ(pos);
2066         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2067         esdTrack->GetESDpid(pid);
2068
2069         vV0FromCascade->AddDaughter(aodTrack =
2070                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2071                                            esdTrack->GetLabel(),
2072                                            p,
2073                                            kTRUE,
2074                                            pos,
2075                                            kFALSE,
2076                                            covTr,
2077                                            (Short_t)esdTrack->Charge(),
2078                                            esdTrack->GetITSClusterMap(),
2079                                            pid,
2080                                            vV0FromCascade,
2081                                            kTRUE,  // check if this is right
2082                                            kFALSE, // check if this is right
2083                                            AliAODTrack::kSecondary)
2084                 );
2085         aodTrack->ConvertAliPIDtoAODPID();
2086       }
2087       else {
2088         cerr << "Error: event " << iEvent << " cascade " << nCascade
2089              << " track " << negFromV0 << " has already been used!" << endl;
2090       }
2091
2092       // Add the bachelor track from the cascade
2093
2094       Int_t bachelor = cascade->GetBindex();
2095
2096       if(!usedTrack[bachelor]) {
2097
2098         usedTrack[bachelor] = kTRUE;
2099
2100         AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2101         esdTrack->GetPxPyPz(p);
2102         esdTrack->GetXYZ(pos);
2103         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2104         esdTrack->GetESDpid(pid);
2105
2106         vcascade->AddDaughter(aodTrack =
2107                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2108                                            esdTrack->GetLabel(),
2109                                            p,
2110                                            kTRUE,
2111                                            pos,
2112                                            kFALSE,
2113                                            covTr,
2114                                            (Short_t)esdTrack->Charge(),
2115                                            esdTrack->GetITSClusterMap(),
2116                                            pid,
2117                                            vcascade,
2118                                            kTRUE,  // check if this is right
2119                                            kFALSE, // check if this is right
2120                                            AliAODTrack::kSecondary)
2121                 );
2122         aodTrack->ConvertAliPIDtoAODPID();
2123      }
2124       else {
2125         cerr << "Error: event " << iEvent << " cascade " << nCascade
2126              << " track " << bachelor << " has already been used!" << endl;
2127       }
2128
2129       // Add the primary track of the cascade (if any)
2130
2131     } // end of the loop on cascades
2132
2133     // V0s
2134
2135     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2136
2137       if (usedV0[nV0]) continue; // skip if aready added to the AOD
2138
2139       AliESDv0 *v0 = esd->GetV0(nV0);
2140
2141       v0->GetXYZ(pos[0], pos[1], pos[2]);
2142       v0->GetPosCov(covVtx);
2143
2144       AliAODVertex * vV0 =
2145         new(vertices[jVertices++]) AliAODVertex(pos,
2146                                                 covVtx,
2147                                                 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2148                                                 primary,
2149                                                 nV0,
2150                                                 AliAODVertex::kV0);
2151       primary->AddDaughter(vV0);
2152
2153       Int_t posFromV0 = v0->GetPindex();
2154       Int_t negFromV0 = v0->GetNindex();
2155
2156       // Add the positive tracks from the V0
2157
2158       if (!usedTrack[posFromV0]) {
2159
2160         usedTrack[posFromV0] = kTRUE;
2161
2162         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2163         esdTrack->GetPxPyPz(p);
2164         esdTrack->GetXYZ(pos);
2165         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2166         esdTrack->GetESDpid(pid);
2167
2168         vV0->AddDaughter(aodTrack =
2169                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2170                                            esdTrack->GetLabel(),
2171                                            p,
2172                                            kTRUE,
2173                                            pos,
2174                                            kFALSE,
2175                                            covTr,
2176                                            (Short_t)esdTrack->Charge(),
2177                                            esdTrack->GetITSClusterMap(),
2178                                            pid,
2179                                            vV0,
2180                                            kTRUE,  // check if this is right
2181                                            kFALSE, // check if this is right
2182                                            AliAODTrack::kSecondary)
2183                 );
2184         aodTrack->ConvertAliPIDtoAODPID();
2185       }
2186       else {
2187         cerr << "Error: event " << iEvent << " V0 " << nV0
2188              << " track " << posFromV0 << " has already been used!" << endl;
2189       }
2190
2191       // Add the negative tracks from the V0
2192
2193       if (!usedTrack[negFromV0]) {
2194
2195         usedTrack[negFromV0] = kTRUE;
2196
2197         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2198         esdTrack->GetPxPyPz(p);
2199         esdTrack->GetXYZ(pos);
2200         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2201         esdTrack->GetESDpid(pid);
2202
2203         vV0->AddDaughter(aodTrack =
2204                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2205                                            esdTrack->GetLabel(),
2206                                            p,
2207                                            kTRUE,
2208                                            pos,
2209                                            kFALSE,
2210                                            covTr,
2211                                            (Short_t)esdTrack->Charge(),
2212                                            esdTrack->GetITSClusterMap(),
2213                                            pid,
2214                                            vV0,
2215                                            kTRUE,  // check if this is right
2216                                            kFALSE, // check if this is right
2217                                            AliAODTrack::kSecondary)
2218                 );
2219         aodTrack->ConvertAliPIDtoAODPID();
2220       }
2221       else {
2222         cerr << "Error: event " << iEvent << " V0 " << nV0
2223              << " track " << negFromV0 << " has already been used!" << endl;
2224       }
2225
2226     } // end of the loop on V0s
2227
2228     // Kinks: it is a big mess the access to the information in the kinks
2229     // The loop is on the tracks in order to find the mother and daugther of each kink
2230
2231
2232     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2233
2234
2235       AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2236
2237       Int_t ikink = esdTrack->GetKinkIndex(0);
2238
2239       if (ikink) {
2240         // Negative kink index: mother, positive: daughter
2241
2242         // Search for the second track of the kink
2243
2244         for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2245
2246           AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2247
2248           Int_t jkink = esdTrack1->GetKinkIndex(0);
2249
2250           if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2251
2252             // The two tracks are from the same kink
2253
2254             if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2255
2256             Int_t imother = -1;
2257             Int_t idaughter = -1;
2258
2259             if (ikink<0 && jkink>0) {
2260
2261               imother = iTrack;
2262               idaughter = jTrack;
2263             }
2264             else if (ikink>0 && jkink<0) {
2265
2266               imother = jTrack;
2267               idaughter = iTrack;
2268             }
2269             else {
2270               cerr << "Error: Wrong combination of kink indexes: "
2271               << ikink << " " << jkink << endl;
2272               continue;
2273             }
2274
2275             // Add the mother track
2276
2277             AliAODTrack * mother = NULL;
2278
2279             if (!usedTrack[imother]) {
2280
2281               usedTrack[imother] = kTRUE;
2282
2283               AliESDtrack *esdTrack = esd->GetTrack(imother);
2284               esdTrack->GetPxPyPz(p);
2285               esdTrack->GetXYZ(pos);
2286               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2287               esdTrack->GetESDpid(pid);
2288
2289               mother =
2290                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2291                                            esdTrack->GetLabel(),
2292                                            p,
2293                                            kTRUE,
2294                                            pos,
2295                                            kFALSE,
2296                                            covTr,
2297                                            (Short_t)esdTrack->Charge(),
2298                                            esdTrack->GetITSClusterMap(),
2299                                            pid,
2300                                            primary,
2301                                            kTRUE, // check if this is right
2302                                            kTRUE, // check if this is right
2303                                            AliAODTrack::kPrimary);
2304               primary->AddDaughter(mother);
2305               mother->ConvertAliPIDtoAODPID();
2306             }
2307             else {
2308               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2309               << " track " << imother << " has already been used!" << endl;
2310             }
2311
2312             // Add the kink vertex
2313             AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2314
2315             AliAODVertex * vkink =
2316             new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2317                                                     NULL,
2318                                                     0.,
2319                                                     mother,
2320                                                     esdTrack->GetID(), // This is the track ID of the mother's track!
2321                                                     AliAODVertex::kKink);
2322             // Add the daughter track
2323
2324             AliAODTrack * daughter = NULL;
2325
2326             if (!usedTrack[idaughter]) {
2327
2328               usedTrack[idaughter] = kTRUE;
2329
2330               AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2331               esdTrack->GetPxPyPz(p);
2332               esdTrack->GetXYZ(pos);
2333               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2334               esdTrack->GetESDpid(pid);
2335
2336               daughter =
2337                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2338                                            esdTrack->GetLabel(),
2339                                            p,
2340                                            kTRUE,
2341                                            pos,
2342                                            kFALSE,
2343                                            covTr,
2344                                            (Short_t)esdTrack->Charge(),
2345                                            esdTrack->GetITSClusterMap(),
2346                                            pid,
2347                                            vkink,
2348                                            kTRUE, // check if this is right
2349                                            kTRUE, // check if this is right
2350                                            AliAODTrack::kPrimary);
2351               vkink->AddDaughter(daughter);
2352               daughter->ConvertAliPIDtoAODPID();
2353             }
2354             else {
2355               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2356               << " track " << idaughter << " has already been used!" << endl;
2357             }
2358
2359
2360           }
2361         }
2362
2363       }
2364
2365     }
2366
2367
2368     // Tracks (primary and orphan)
2369
2370     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2371
2372
2373       if (usedTrack[nTrack]) continue;
2374
2375       AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2376       esdTrack->GetPxPyPz(p);
2377       esdTrack->GetXYZ(pos);
2378       esdTrack->GetCovarianceXYZPxPyPz(covTr);
2379       esdTrack->GetESDpid(pid);
2380
2381       Float_t impactXY, impactZ;
2382
2383       esdTrack->GetImpactParameters(impactXY,impactZ);
2384
2385       if (impactXY<3) {
2386         // track inside the beam pipe
2387
2388         primary->AddDaughter(aodTrack =
2389             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2390                                          esdTrack->GetLabel(),
2391                                          p,
2392                                          kTRUE,
2393                                          pos,
2394                                          kFALSE,
2395                                          covTr,
2396                                          (Short_t)esdTrack->Charge(),
2397                                          esdTrack->GetITSClusterMap(),
2398                                          pid,
2399                                          primary,
2400                                          kTRUE, // check if this is right
2401                                          kTRUE, // check if this is right
2402                                          AliAODTrack::kPrimary)
2403             );
2404         aodTrack->ConvertAliPIDtoAODPID();
2405       }
2406       else {
2407         // outside the beam pipe: orphan track
2408             aodTrack =
2409             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2410                                          esdTrack->GetLabel(),
2411                                          p,
2412                                          kTRUE,
2413                                          pos,
2414                                          kFALSE,
2415                                          covTr,
2416                                          (Short_t)esdTrack->Charge(),
2417                                          esdTrack->GetITSClusterMap(),
2418                                          pid,
2419                                          NULL,
2420                                          kFALSE, // check if this is right
2421                                          kFALSE, // check if this is right
2422                                          AliAODTrack::kOrphan);
2423             aodTrack->ConvertAliPIDtoAODPID();
2424       }
2425     } // end of loop on tracks
2426
2427     // muon tracks
2428     Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2429     for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2430
2431       AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2432       p[0] = esdMuTrack->Px();
2433       p[1] = esdMuTrack->Py();
2434       p[2] = esdMuTrack->Pz();
2435       pos[0] = primary->GetX();
2436       pos[1] = primary->GetY();
2437       pos[2] = primary->GetZ();
2438
2439       // has to be changed once the muon pid is provided by the ESD
2440       for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2441
2442       primary->AddDaughter(aodTrack =
2443           new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2444                                              0, // no label provided
2445                                              p,
2446                                              kTRUE,
2447                                              pos,
2448                                              kFALSE,
2449                                              NULL, // no covariance matrix provided
2450                                              esdMuTrack->Charge(),
2451                                              0, // ITSClusterMap is set below
2452                                              pid,
2453                                              primary,
2454                                              kFALSE,  // muon tracks are not used to fit the primary vtx
2455                                              kFALSE,  // not used for vertex fit
2456                                              AliAODTrack::kPrimary)
2457           );
2458
2459       aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
2460       Int_t track2Trigger = esdMuTrack->GetMatchTrigger();
2461       aodTrack->SetMatchTrigger(track2Trigger);
2462       if (track2Trigger)
2463         aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
2464       else
2465         aodTrack->SetChi2MatchTrigger(0.);
2466     }
2467
2468     // Access to the AOD container of clusters
2469     TClonesArray &clusters = *(aod->GetClusters());
2470     Int_t jClusters=0;
2471
2472     // Calo Clusters
2473     Int_t nClusters    = esd->GetNumberOfCaloClusters();
2474
2475     for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2476
2477       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2478
2479       Int_t id = cluster->GetID();
2480       Int_t label = -1;
2481       Float_t energy = cluster->E();
2482       cluster->GetPosition(posF);
2483       AliAODVertex *prodVertex = primary;
2484       AliAODTrack *primTrack = NULL;
2485       Char_t ttype=AliAODCluster::kUndef;
2486
2487       if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2488       else if (cluster->IsEMCAL()) {
2489
2490         if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2491           ttype = AliAODCluster::kEMCALPseudoCluster;
2492         else
2493           ttype = AliAODCluster::kEMCALClusterv1;
2494
2495       }
2496
2497       new(clusters[jClusters++]) AliAODCluster(id,
2498                                                label,
2499                                                energy,
2500                                                pos,
2501                                                NULL, // no covariance matrix provided
2502                                                NULL, // no pid for clusters provided
2503                                                prodVertex,
2504                                                primTrack,
2505                                                ttype);
2506
2507     } // end of loop on calo clusters
2508
2509     // tracklets
2510     const AliMultiplicity *mult = esd->GetMultiplicity();
2511     if (mult) {
2512       if (mult->GetNumberOfTracklets()>0) {
2513         aod->GetTracklets()->CreateContainer(mult->GetNumberOfTracklets());
2514
2515         for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
2516           aod->GetTracklets()->SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n));
2517         }
2518       }
2519     } else {
2520       Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
2521     }
2522
2523     delete [] usedTrack;
2524     delete [] usedV0;
2525     delete [] usedKink;
2526
2527     // fill the tree for this event
2528     aodTree->Fill();
2529   } // end of event loop
2530
2531   aodTree->GetUserInfo()->Add(aod);
2532
2533   // close ESD file
2534   esdFile->Close();
2535
2536   // write the tree to the specified file
2537   aodFile = aodTree->GetCurrentFile();
2538   aodFile->cd();
2539   aodTree->Write();
2540
2541   return;
2542 }
2543
2544 void AliReconstruction::WriteAlignmentData(AliESDEvent* esd)
2545 {
2546   // Write space-points which are then used in the alignment procedures
2547   // For the moment only ITS, TRD and TPC
2548
2549   // Load TOF clusters
2550   if (fTracker[3]){
2551     fLoader[3]->LoadRecPoints("read");
2552     TTree* tree = fLoader[3]->TreeR();
2553     if (!tree) {
2554       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2555       return;
2556     }
2557     fTracker[3]->LoadClusters(tree);
2558   }
2559   Int_t ntracks = esd->GetNumberOfTracks();
2560   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2561     {
2562       AliESDtrack *track = esd->GetTrack(itrack);
2563       Int_t nsp = 0;
2564       Int_t idx[200];
2565       for (Int_t iDet = 3; iDet >= 0; iDet--)
2566         nsp += track->GetNcls(iDet);
2567       if (nsp) {
2568         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2569         track->SetTrackPointArray(sp);
2570         Int_t isptrack = 0;
2571         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2572           AliTracker *tracker = fTracker[iDet];
2573           if (!tracker) continue;
2574           Int_t nspdet = track->GetNcls(iDet);
2575           if (nspdet <= 0) continue;
2576           track->GetClusters(iDet,idx);
2577           AliTrackPoint p;
2578           Int_t isp = 0;
2579           Int_t isp2 = 0;
2580           while (isp < nspdet) {
2581             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2582             const Int_t kNTPCmax = 159;
2583             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2584             if (!isvalid) continue;
2585             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2586           }
2587         }
2588       }
2589     }
2590   if (fTracker[3]){
2591     fTracker[3]->UnloadClusters();
2592     fLoader[3]->UnloadRecPoints();
2593   }
2594 }
2595
2596 //_____________________________________________________________________________
2597 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESDEvent* esd)
2598 {
2599   // The method reads the raw-data error log
2600   // accumulated within the rawReader.
2601   // It extracts the raw-data errors related to
2602   // the current event and stores them into
2603   // a TClonesArray inside the esd object.
2604
2605   if (!fRawReader) return;
2606
2607   for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2608
2609     AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2610     if (!log) continue;
2611     if (iEvent != log->GetEventNumber()) continue;
2612
2613     esd->AddRawDataErrorLog(log);
2614   }
2615
2616 }
2617
2618 TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2619   // Dump a file content into a char in TNamed
2620   ifstream in;
2621   in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2622   Int_t kBytes = (Int_t)in.tellg();
2623   printf("Size: %d \n",kBytes);
2624   TNamed *fn = 0;
2625   if(in.good()){
2626     char* memblock = new char [kBytes];
2627     in.seekg (0, ios::beg);
2628     in.read (memblock, kBytes);
2629     in.close();
2630     TString fData(memblock,kBytes);
2631     fn = new TNamed(fName,fData);
2632     printf("fData Size: %d \n",fData.Sizeof());
2633     printf("fName Size: %d \n",fName.Sizeof());
2634     printf("fn    Size: %d \n",fn->Sizeof());
2635     delete[] memblock;
2636   }
2637   else{
2638     AliInfo(Form("Could not Open %s\n",fPath.Data()));
2639   }
2640
2641   return fn;
2642 }
2643
2644 void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2645   // This is not really needed in AliReconstruction at the moment
2646   // but can serve as a template
2647
2648   TList *fList = fTree->GetUserInfo();
2649   TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2650   printf("fn Size: %d \n",fn->Sizeof());
2651
2652   TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2653   const char* cdata = fn->GetTitle();
2654   printf("fTmp Size %d\n",fTmp.Sizeof());
2655
2656   int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2657   printf("calculated size %d\n",size);
2658   ofstream out(fName.Data(),ios::out | ios::binary);
2659   out.write(cdata,size);
2660   out.close();
2661
2662 }
2663
2664 //_____________________________________________________________________________
2665 AliQualAssDataMaker * AliReconstruction::GetQualAssDataMaker(Int_t iDet)
2666 {
2667 // get the quality assurance data maker object and the loader for a detector
2668
2669   if (fQualAssDataMaker[iDet])
2670     return fQualAssDataMaker[iDet];
2671
2672   // load the QA data maker object
2673   TPluginManager* pluginManager = gROOT->GetPluginManager();
2674   TString detName = fgkDetectorName[iDet];
2675   TString qadmName = "Ali" + detName + "QualAssDataMaker";
2676   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT"))
2677     return NULL;
2678
2679   AliQualAssDataMaker * qadm = NULL;
2680   // first check if a plugin is defined for the quality assurance data maker
2681   TPluginHandler* pluginHandler = pluginManager->FindHandler("AliQualAssDataMaker", detName);
2682   // if not, add a plugin for it
2683   if (!pluginHandler) {
2684     AliDebug(1, Form("defining plugin for %s", qadmName.Data()));
2685     TString libs = gSystem->GetLibraries();
2686     if (libs.Contains("lib" + detName + "base.so") ||
2687         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
2688       pluginManager->AddHandler("AliQualAssDataMaker", detName,
2689                                 qadmName, detName + "qadm", qadmName + "()");
2690     } else {
2691       pluginManager->AddHandler("AliQualAssDataMaker", detName,
2692                                 qadmName, detName, qadmName + "()");
2693     }
2694     pluginHandler = pluginManager->FindHandler("AliQualAssDataMaker", detName);
2695   }
2696   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
2697     qadm = (AliQualAssDataMaker *) pluginHandler->ExecPlugin(0);
2698   }
2699   if (qadm) {
2700     AliInfo(Form("Initializing quality assurance data maker for %s", fgkDetectorName[iDet]));
2701     qadm->Init(AliQualAss::kRECPOINTS);
2702     qadm->Init(AliQualAss::kESDS) ;
2703     fQualAssDataMaker[iDet] = qadm;
2704   }
2705
2706   return qadm;
2707 }
2708
2709 //_____________________________________________________________________________
2710 Bool_t AliReconstruction::RunQualAss(const char* detectors, AliESDEvent *& esd)
2711 {
2712   // run the Quality Assurance data producer
2713
2714   AliCodeTimerAuto("")
2715   TString detStr = detectors;
2716   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
2717    if (!IsSelected(fgkDetectorName[iDet], detStr))
2718      continue;
2719    AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
2720    if (!qadm)
2721      continue;
2722    AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2723    AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2724
2725    qadm->Exec(AliQualAss::kESDS, esd) ;
2726
2727    AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2728  }
2729  if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
2730    AliError(Form("the following detectors were not found: %s",
2731                  detStr.Data()));
2732    if (fStopOnError)
2733      return kFALSE;
2734  }
2735
2736  return kTRUE;
2737 }