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