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