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