STEER/AliReconstruction.cxx

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