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->GetAliRun()->GetRunNumber());
 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->LoadKinematics();
1610
1611   } else {               // galice.root does not exist
1612     if (!fRawReader) {
1613       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1614       CleanUp();
1615       return kFALSE;
1616     }
1617     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1618                                     AliConfig::GetDefaultEventFolderName(),
1619                                     "recreate");
1620     if (!fRunLoader) {
1621       AliError(Form("could not create run loader in file %s",
1622                     fGAliceFileName.Data()));
1623       CleanUp();
1624       return kFALSE;
1625     }
1626     fRunLoader->MakeTree("E");
1627     Int_t iEvent = 0;
1628     while (fRawReader->NextEvent()) {
1629       fRunLoader->SetEventNumber(iEvent);
1630       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1631                                      iEvent, iEvent);
1632       fRunLoader->MakeTree("H");
1633       fRunLoader->TreeE()->Fill();
1634       iEvent++;
1635     }
1636     fRawReader->RewindEvents();
1637     if (fNumberOfEventsPerFile > 0)
1638       fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1639     else
1640       fRunLoader->SetNumberOfEventsPerFile(iEvent);
1641     fRunLoader->WriteHeader("OVERWRITE");
1642     fRunLoader->CdGAFile();
1643     fRunLoader->Write(0, TObject::kOverwrite);
1644 //    AliTracker::SetFieldMap(???);
1645   }
1646
1647   return kTRUE;
1648 }
1649
1650 //_____________________________________________________________________________
1651 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1652 {
1653 // get the reconstructor object and the loader for a detector
1654
1655   if (fReconstructor[iDet]) return fReconstructor[iDet];
1656
1657   // load the reconstructor object
1658   TPluginManager* pluginManager = gROOT->GetPluginManager();
1659   TString detName = fgkDetectorName[iDet];
1660   TString recName = "Ali" + detName + "Reconstructor";
1661   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1662
1663   AliReconstructor* reconstructor = NULL;
1664   // first check if a plugin is defined for the reconstructor
1665   TPluginHandler* pluginHandler =
1666     pluginManager->FindHandler("AliReconstructor", detName);
1667   // if not, add a plugin for it
1668   if (!pluginHandler) {
1669     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1670     TString libs = gSystem->GetLibraries();
1671     if (libs.Contains("lib" + detName + "base.so") ||
1672         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1673       pluginManager->AddHandler("AliReconstructor", detName,
1674                                 recName, detName + "rec", recName + "()");
1675     } else {
1676       pluginManager->AddHandler("AliReconstructor", detName,
1677                                 recName, detName, recName + "()");
1678     }
1679     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1680   }
1681   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1682     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1683   }
1684   if (reconstructor) {
1685     TObject* obj = fOptions.FindObject(detName.Data());
1686     if (obj) reconstructor->SetOption(obj->GetTitle());
1687     reconstructor->Init();
1688     fReconstructor[iDet] = reconstructor;
1689   }
1690
1691   // get or create the loader
1692   if (detName != "HLT") {
1693     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1694     if (!fLoader[iDet]) {
1695       AliConfig::Instance()
1696         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1697                                 detName, detName);
1698       // first check if a plugin is defined for the loader
1699       pluginHandler =
1700         pluginManager->FindHandler("AliLoader", detName);
1701       // if not, add a plugin for it
1702       if (!pluginHandler) {
1703         TString loaderName = "Ali" + detName + "Loader";
1704         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1705         pluginManager->AddHandler("AliLoader", detName,
1706                                   loaderName, detName + "base",
1707                                   loaderName + "(const char*, TFolder*)");
1708         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1709       }
1710       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1711         fLoader[iDet] =
1712           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1713                                                  fRunLoader->GetEventFolder());
1714       }
1715       if (!fLoader[iDet]) {   // use default loader
1716         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1717       }
1718       if (!fLoader[iDet]) {
1719         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1720         if (fStopOnError) return NULL;
1721       } else {
1722         fRunLoader->AddLoader(fLoader[iDet]);
1723         fRunLoader->CdGAFile();
1724         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1725         fRunLoader->Write(0, TObject::kOverwrite);
1726       }
1727     }
1728   }
1729
1730   return reconstructor;
1731 }
1732
1733 //_____________________________________________________________________________
1734 Bool_t AliReconstruction::CreateVertexer()
1735 {
1736 // create the vertexer
1737
1738   fVertexer = NULL;
1739   AliReconstructor* itsReconstructor = GetReconstructor(0);
1740   if (itsReconstructor) {
1741     fVertexer = itsReconstructor->CreateVertexer();
1742   }
1743   if (!fVertexer) {
1744     AliWarning("couldn't create a vertexer for ITS");
1745     if (fStopOnError) return kFALSE;
1746   }
1747
1748   return kTRUE;
1749 }
1750
1751 //_____________________________________________________________________________
1752 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1753 {
1754 // create the trackers
1755
1756   TString detStr = detectors;
1757   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1758     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1759     AliReconstructor* reconstructor = GetReconstructor(iDet);
1760     if (!reconstructor) continue;
1761     TString detName = fgkDetectorName[iDet];
1762     if (detName == "HLT") {
1763       fRunHLTTracking = kTRUE;
1764       continue;
1765     }
1766     if (detName == "MUON") {
1767       fRunMuonTracking = kTRUE;
1768       continue;
1769     }
1770
1771
1772     fTracker[iDet] = reconstructor->CreateTracker();
1773     if (!fTracker[iDet] && (iDet < 7)) {
1774       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1775       if (fStopOnError) return kFALSE;
1776     }
1777     AliSysInfo::AddStamp(Form("LTracker%s",fgkDetectorName[iDet]), iDet,0);
1778   }
1779
1780   return kTRUE;
1781 }
1782
1783 //_____________________________________________________________________________
1784 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1785 {
1786 // delete trackers and the run loader and close and delete the file
1787
1788   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1789     delete fReconstructor[iDet];
1790     fReconstructor[iDet] = NULL;
1791     fLoader[iDet] = NULL;
1792     delete fTracker[iDet];
1793     fTracker[iDet] = NULL;
1794     delete fQADataMaker[iDet];
1795     fQADataMaker[iDet] = NULL;
1796   }
1797   delete fVertexer;
1798   fVertexer = NULL;
1799   delete fDiamondProfile;
1800   fDiamondProfile = NULL;
1801
1802   delete fRunLoader;
1803   fRunLoader = NULL;
1804   delete fRawReader;
1805   fRawReader = NULL;
1806
1807   if (file) {
1808     file->Close();
1809     delete file;
1810   }
1811
1812   if (fileOld) {
1813     fileOld->Close();
1814     delete fileOld;
1815     gSystem->Unlink("AliESDs.old.root");
1816   }
1817 }
1818
1819 //_____________________________________________________________________________
1820
1821 Bool_t AliReconstruction::ReadESD(AliESDEvent*& esd, const char* recStep) const
1822 {
1823 // read the ESD event from a file
1824
1825   if (!esd) return kFALSE;
1826   char fileName[256];
1827   sprintf(fileName, "ESD_%d.%d_%s.root",
1828           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1829   if (gSystem->AccessPathName(fileName)) return kFALSE;
1830
1831   AliInfo(Form("reading ESD from file %s", fileName));
1832   AliDebug(1, Form("reading ESD from file %s", fileName));
1833   TFile* file = TFile::Open(fileName);
1834   if (!file || !file->IsOpen()) {
1835     AliError(Form("opening %s failed", fileName));
1836     delete file;
1837     return kFALSE;
1838   }
1839
1840   gROOT->cd();
1841   delete esd;
1842   esd = (AliESDEvent*) file->Get("ESD");
1843   file->Close();
1844   delete file;
1845   return kTRUE;
1846
1847 }
1848
1849
1850
1851 //_____________________________________________________________________________
1852 void AliReconstruction::WriteESD(AliESDEvent* esd, const char* recStep) const
1853 {
1854 // write the ESD event to a file
1855
1856   if (!esd) return;
1857   char fileName[256];
1858   sprintf(fileName, "ESD_%d.%d_%s.root",
1859           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1860
1861   AliDebug(1, Form("writing ESD to file %s", fileName));
1862   TFile* file = TFile::Open(fileName, "recreate");
1863   if (!file || !file->IsOpen()) {
1864     AliError(Form("opening %s failed", fileName));
1865   } else {
1866     esd->Write("ESD");
1867     file->Close();
1868   }
1869   delete file;
1870 }
1871
1872
1873
1874
1875
1876 //_____________________________________________________________________________
1877 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
1878 {
1879   // write all files from the given esd file to an aod file
1880
1881   // create an AliAOD object
1882   AliAODEvent *aod = new AliAODEvent();
1883   aod->CreateStdContent();
1884
1885   // go to the file
1886   aodFile->cd();
1887
1888   // create the tree
1889   TTree *aodTree = new TTree("aodTree", "AliAOD tree");
1890   aodTree->Branch(aod->GetList());
1891
1892   // connect to ESD
1893   TTree *t = (TTree*) esdFile->Get("esdTree");
1894   AliESDEvent *esd = new AliESDEvent();
1895   esd->ReadFromTree(t);
1896
1897   Int_t nEvents = t->GetEntries();
1898
1899   // set arrays and pointers
1900   Float_t posF[3];
1901   Double_t pos[3];
1902   Double_t p[3];
1903   Double_t p_pos[3];
1904   Double_t p_neg[3];
1905   Double_t covVtx[6];
1906   Double_t covTr[21];
1907   Double_t pid[10];
1908
1909   // loop over events and fill them
1910   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
1911     //cout << "event: " << iEvent << endl;
1912     t->GetEntry(iEvent);
1913
1914     // Multiplicity information needed by the header (to be revised!)
1915     Int_t nTracks   = esd->GetNumberOfTracks();
1916     Int_t nPosTracks = 0;
1917     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
1918       if (esd->GetTrack(iTrack)->Charge()> 0) nPosTracks++;
1919
1920     // Access the header
1921     AliAODHeader *header = aod->GetHeader();
1922
1923     // fill the header
1924     header->SetRunNumber       (esd->GetRunNumber()       );
1925     header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
1926     header->SetOrbitNumber     (esd->GetOrbitNumber()     );
1927     header->SetPeriodNumber    (esd->GetPeriodNumber()    );
1928     header->SetTriggerMask     (esd->GetTriggerMask()     );
1929     header->SetTriggerCluster  (esd->GetTriggerCluster()  );
1930     header->SetEventType       (esd->GetEventType()       );
1931     header->SetMagneticField   (esd->GetMagneticField()   );
1932     header->SetZDCN1Energy     (esd->GetZDCN1Energy()     );
1933     header->SetZDCP1Energy     (esd->GetZDCP1Energy()     );
1934     header->SetZDCN2Energy     (esd->GetZDCN2Energy()     );
1935     header->SetZDCP2Energy     (esd->GetZDCP2Energy()     );
1936     header->SetZDCEMEnergy     (esd->GetZDCEMEnergy()     );
1937     header->SetRefMultiplicity   (nTracks);
1938     header->SetRefMultiplicityPos(nPosTracks);
1939     header->SetRefMultiplicityNeg(nTracks - nPosTracks);
1940     header->SetMuonMagFieldScale(-999.); // FIXME
1941     header->SetCentrality(-999.);        // FIXME
1942
1943     Int_t nV0s      = esd->GetNumberOfV0s();
1944     Int_t nCascades = esd->GetNumberOfCascades();
1945     Int_t nKinks    = esd->GetNumberOfKinks();
1946     Int_t nVertices = nV0s + nCascades + nKinks + 1 /* = prim. vtx*/;
1947     Int_t nJets     = 0;
1948     Int_t nCaloClus = esd->GetNumberOfCaloClusters();
1949     Int_t nFmdClus  = 0;
1950     Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
1951
1952     aod->ResetStd(nTracks, nVertices, nV0s+nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
1953
1954     // Array to take into account the tracks already added to the AOD
1955     Bool_t * usedTrack = NULL;
1956     if (nTracks>0) {
1957       usedTrack = new Bool_t[nTracks];
1958       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
1959     }
1960     // Array to take into account the V0s already added to the AOD
1961     Bool_t * usedV0 = NULL;
1962     if (nV0s>0) {
1963       usedV0 = new Bool_t[nV0s];
1964       for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
1965     }
1966     // Array to take into account the kinks already added to the AOD
1967     Bool_t * usedKink = NULL;
1968     if (nKinks>0) {
1969       usedKink = new Bool_t[nKinks];
1970       for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
1971     }
1972
1973     // Access to the AOD container of vertices
1974     TClonesArray &vertices = *(aod->GetVertices());
1975     Int_t jVertices=0;
1976
1977     // Access to the AOD container of tracks
1978     TClonesArray &tracks = *(aod->GetTracks());
1979     Int_t jTracks=0;
1980
1981     // Access to the AOD container of V0s
1982     TClonesArray &V0s = *(aod->GetV0s());
1983     Int_t jV0s=0;
1984
1985     // Add primary vertex. The primary tracks will be defined
1986     // after the loops on the composite objects (V0, cascades, kinks)
1987     const AliESDVertex *vtx = esd->GetPrimaryVertex();
1988
1989     vtx->GetXYZ(pos); // position
1990     vtx->GetCovMatrix(covVtx); //covariance matrix
1991
1992     AliAODVertex * primary = new(vertices[jVertices++])
1993       AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1994
1995
1996     AliAODTrack *aodTrack = 0x0;
1997
1998     // Create vertices starting from the most complex objects
1999
2000     // Cascades
2001     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2002       AliESDcascade *cascade = esd->GetCascade(nCascade);
2003
2004       cascade->GetXYZ(pos[0], pos[1], pos[2]);
2005       cascade->GetPosCovXi(covVtx);
2006
2007       // Add the cascade vertex
2008       AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2009                                                                         covVtx,
2010                                                                         cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2011                                                                         primary,
2012                                                                         nCascade,
2013                                                                         AliAODVertex::kCascade);
2014
2015       primary->AddDaughter(vcascade); // the cascade 'particle' (represented by a vertex) is added as a daughter to the primary vertex
2016
2017       // Add the V0 from the cascade. The ESD class have to be optimized...
2018       // Now we have to search for the corresponding V0 in the list of V0s
2019       // using the indeces of the positive and negative tracks
2020
2021       Int_t posFromV0 = cascade->GetPindex();
2022       Int_t negFromV0 = cascade->GetNindex();
2023
2024
2025       AliESDv0 * v0 = 0x0;
2026       Int_t indV0 = -1;
2027
2028       for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2029
2030         v0 = esd->GetV0(iV0);
2031         Int_t posV0 = v0->GetPindex();
2032         Int_t negV0 = v0->GetNindex();
2033
2034         if (posV0==posFromV0 && negV0==negFromV0) {
2035           indV0 = iV0;
2036           break;
2037         }
2038       }
2039
2040       AliAODVertex * vV0FromCascade = 0x0;
2041
2042       if (indV0>-1 && !usedV0[indV0]) {
2043
2044         // the V0 exists in the array of V0s and is not used
2045
2046         usedV0[indV0] = kTRUE;
2047
2048         v0->GetXYZ(pos[0], pos[1], pos[2]);
2049         v0->GetPosCov(covVtx);
2050
2051         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2052                                                                  covVtx,
2053                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2054                                                                  vcascade,
2055                                                                  indV0,
2056                                                                  AliAODVertex::kV0);
2057       } else {
2058
2059         // the V0 doesn't exist in the array of V0s or was used
2060         cerr << "Error: event " << iEvent << " cascade " << nCascade
2061              << " The V0 " << indV0
2062              << " doesn't exist in the array of V0s or was used!" << endl;
2063
2064         cascade->GetXYZ(pos[0], pos[1], pos[2]);
2065         cascade->GetPosCov(covVtx);
2066
2067         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2068                                                                  covVtx,
2069                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2070                                                                  vcascade,
2071                                                                  indV0,
2072                                                                  AliAODVertex::kV0);
2073         vcascade->AddDaughter(vV0FromCascade);
2074
2075       }
2076
2077       // Add the positive tracks from the V0
2078
2079       if (! usedTrack[posFromV0]) {
2080
2081         usedTrack[posFromV0] = kTRUE;
2082
2083         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2084         esdTrack->GetPxPyPz(p_pos);
2085         esdTrack->GetXYZ(pos);
2086         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2087         esdTrack->GetESDpid(pid);
2088
2089         vV0FromCascade->AddDaughter(aodTrack =
2090                                     new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2091                                            esdTrack->GetLabel(),
2092                                            p_pos,
2093                                            kTRUE,
2094                                            pos,
2095                                            kFALSE,
2096                                            covTr,
2097                                            (Short_t)esdTrack->Charge(),
2098                                            esdTrack->GetITSClusterMap(),
2099                                            pid,
2100                                            vV0FromCascade,
2101                                            kTRUE,  // check if this is right
2102                                            kFALSE, // check if this is right
2103                                            AliAODTrack::kSecondary)
2104                 );
2105         aodTrack->ConvertAliPIDtoAODPID();
2106       }
2107       else {
2108         cerr << "Error: event " << iEvent << " cascade " << nCascade
2109              << " track " << posFromV0 << " has already been used!" << endl;
2110       }
2111
2112       // Add the negative tracks from the V0
2113
2114       if (!usedTrack[negFromV0]) {
2115
2116         usedTrack[negFromV0] = kTRUE;
2117
2118         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2119         esdTrack->GetPxPyPz(p_neg);
2120         esdTrack->GetXYZ(pos);
2121         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2122         esdTrack->GetESDpid(pid);
2123
2124         vV0FromCascade->AddDaughter(aodTrack =
2125                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2126                                            esdTrack->GetLabel(),
2127                                            p_neg,
2128                                            kTRUE,
2129                                            pos,
2130                                            kFALSE,
2131                                            covTr,
2132                                            (Short_t)esdTrack->Charge(),
2133                                            esdTrack->GetITSClusterMap(),
2134                                            pid,
2135                                            vV0FromCascade,
2136                                            kTRUE,  // check if this is right
2137                                            kFALSE, // check if this is right
2138                                            AliAODTrack::kSecondary)
2139                 );
2140         aodTrack->ConvertAliPIDtoAODPID();
2141       }
2142       else {
2143         cerr << "Error: event " << iEvent << " cascade " << nCascade
2144              << " track " << negFromV0 << " has already been used!" << endl;
2145       }
2146
2147       // add it to the V0 array as well
2148       Double_t d0[2] = { -999., -99.};
2149       // counting is probably wrong
2150       new(V0s[jV0s++]) AliAODv0(vV0FromCascade, -999., -99., p_pos, p_neg, d0); // to be refined
2151
2152       // Add the bachelor track from the cascade
2153
2154       Int_t bachelor = cascade->GetBindex();
2155
2156       if(!usedTrack[bachelor]) {
2157
2158         usedTrack[bachelor] = kTRUE;
2159
2160         AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2161         esdTrack->GetPxPyPz(p);
2162         esdTrack->GetXYZ(pos);
2163         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2164         esdTrack->GetESDpid(pid);
2165
2166         vcascade->AddDaughter(aodTrack =
2167                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2168                                            esdTrack->GetLabel(),
2169                                            p,
2170                                            kTRUE,
2171                                            pos,
2172                                            kFALSE,
2173                                            covTr,
2174                                            (Short_t)esdTrack->Charge(),
2175                                            esdTrack->GetITSClusterMap(),
2176                                            pid,
2177                                            vcascade,
2178                                            kTRUE,  // check if this is right
2179                                            kFALSE, // check if this is right
2180                                            AliAODTrack::kSecondary)
2181                 );
2182         aodTrack->ConvertAliPIDtoAODPID();
2183      }
2184       else {
2185         cerr << "Error: event " << iEvent << " cascade " << nCascade
2186              << " track " << bachelor << " has already been used!" << endl;
2187       }
2188
2189       // Add the primary track of the cascade (if any)
2190
2191     } // end of the loop on cascades
2192
2193     // V0s
2194
2195     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2196
2197       if (usedV0[nV0]) continue; // skip if aready added to the AOD
2198
2199       AliESDv0 *v0 = esd->GetV0(nV0);
2200
2201       v0->GetXYZ(pos[0], pos[1], pos[2]);
2202       v0->GetPosCov(covVtx);
2203
2204       AliAODVertex * vV0 =
2205         new(vertices[jVertices++]) AliAODVertex(pos,
2206                                                 covVtx,
2207                                                 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2208                                                 primary,
2209                                                 nV0,
2210                                                 AliAODVertex::kV0);
2211       primary->AddDaughter(vV0);
2212
2213       Int_t posFromV0 = v0->GetPindex();
2214       Int_t negFromV0 = v0->GetNindex();
2215
2216       // Add the positive tracks from the V0
2217
2218       if (!usedTrack[posFromV0]) {
2219
2220         usedTrack[posFromV0] = kTRUE;
2221
2222         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2223         esdTrack->GetPxPyPz(p_pos);
2224         esdTrack->GetXYZ(pos);
2225         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2226         esdTrack->GetESDpid(pid);
2227
2228         vV0->AddDaughter(aodTrack =
2229                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2230                                            esdTrack->GetLabel(),
2231                                            p_pos,
2232                                            kTRUE,
2233                                            pos,
2234                                            kFALSE,
2235                                            covTr,
2236                                            (Short_t)esdTrack->Charge(),
2237                                            esdTrack->GetITSClusterMap(),
2238                                            pid,
2239                                            vV0,
2240                                            kTRUE,  // check if this is right
2241                                            kFALSE, // check if this is right
2242                                            AliAODTrack::kSecondary)
2243                 );
2244         aodTrack->ConvertAliPIDtoAODPID();
2245       }
2246       else {
2247         cerr << "Error: event " << iEvent << " V0 " << nV0
2248              << " track " << posFromV0 << " has already been used!" << endl;
2249       }
2250
2251       // Add the negative tracks from the V0
2252
2253       if (!usedTrack[negFromV0]) {
2254
2255         usedTrack[negFromV0] = kTRUE;
2256
2257         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2258         esdTrack->GetPxPyPz(p_neg);
2259         esdTrack->GetXYZ(pos);
2260         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2261         esdTrack->GetESDpid(pid);
2262
2263         vV0->AddDaughter(aodTrack =
2264                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2265                                            esdTrack->GetLabel(),
2266                                            p_neg,
2267                                            kTRUE,
2268                                            pos,
2269                                            kFALSE,
2270                                            covTr,
2271                                            (Short_t)esdTrack->Charge(),
2272                                            esdTrack->GetITSClusterMap(),
2273                                            pid,
2274                                            vV0,
2275                                            kTRUE,  // check if this is right
2276                                            kFALSE, // check if this is right
2277                                            AliAODTrack::kSecondary)
2278                 );
2279         aodTrack->ConvertAliPIDtoAODPID();
2280       }
2281       else {
2282         cerr << "Error: event " << iEvent << " V0 " << nV0
2283              << " track " << negFromV0 << " has already been used!" << endl;
2284       }
2285
2286       // add it to the V0 array as well
2287       Double_t d0[2] = { 999., 99.};
2288       new(V0s[jV0s++]) AliAODv0(vV0, 999., 99., p_pos, p_neg, d0); // to be refined
2289     } // end of the loop on V0s
2290
2291     // Kinks: it is a big mess the access to the information in the kinks
2292     // The loop is on the tracks in order to find the mother and daugther of each kink
2293
2294
2295     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2296
2297       AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2298
2299       Int_t ikink = esdTrack->GetKinkIndex(0);
2300
2301       if (ikink) {
2302         // Negative kink index: mother, positive: daughter
2303
2304         // Search for the second track of the kink
2305
2306         for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2307
2308           AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2309
2310           Int_t jkink = esdTrack1->GetKinkIndex(0);
2311
2312           if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2313
2314             // The two tracks are from the same kink
2315
2316             if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2317
2318             Int_t imother = -1;
2319             Int_t idaughter = -1;
2320
2321             if (ikink<0 && jkink>0) {
2322
2323               imother = iTrack;
2324               idaughter = jTrack;
2325             }
2326             else if (ikink>0 && jkink<0) {
2327
2328               imother = jTrack;
2329               idaughter = iTrack;
2330             }
2331             else {
2332               cerr << "Error: Wrong combination of kink indexes: "
2333               << ikink << " " << jkink << endl;
2334               continue;
2335             }
2336
2337             // Add the mother track
2338
2339             AliAODTrack * mother = NULL;
2340
2341             if (!usedTrack[imother]) {
2342
2343               usedTrack[imother] = kTRUE;
2344
2345               AliESDtrack *esdTrack = esd->GetTrack(imother);
2346               esdTrack->GetPxPyPz(p);
2347               esdTrack->GetXYZ(pos);
2348               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2349               esdTrack->GetESDpid(pid);
2350
2351               mother =
2352                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2353                                            esdTrack->GetLabel(),
2354                                            p,
2355                                            kTRUE,
2356                                            pos,
2357                                            kFALSE,
2358                                            covTr,
2359                                            (Short_t)esdTrack->Charge(),
2360                                            esdTrack->GetITSClusterMap(),
2361                                            pid,
2362                                            primary,
2363                                            kTRUE, // check if this is right
2364                                            kTRUE, // check if this is right
2365                                            AliAODTrack::kPrimary);
2366               primary->AddDaughter(mother);
2367               mother->ConvertAliPIDtoAODPID();
2368             }
2369             else {
2370               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2371               << " track " << imother << " has already been used!" << endl;
2372             }
2373
2374             // Add the kink vertex
2375             AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2376
2377             AliAODVertex * vkink =
2378             new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2379                                                     NULL,
2380                                                     0.,
2381                                                     mother,
2382                                                     esdTrack->GetID(), // This is the track ID of the mother's track!
2383                                                     AliAODVertex::kKink);
2384             // Add the daughter track
2385
2386             AliAODTrack * daughter = NULL;
2387
2388             if (!usedTrack[idaughter]) {
2389
2390               usedTrack[idaughter] = kTRUE;
2391
2392               AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2393               esdTrack->GetPxPyPz(p);
2394               esdTrack->GetXYZ(pos);
2395               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2396               esdTrack->GetESDpid(pid);
2397
2398               daughter =
2399                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2400                                            esdTrack->GetLabel(),
2401                                            p,
2402                                            kTRUE,
2403                                            pos,
2404                                            kFALSE,
2405                                            covTr,
2406                                            (Short_t)esdTrack->Charge(),
2407                                            esdTrack->GetITSClusterMap(),
2408                                            pid,
2409                                            vkink,
2410                                            kTRUE, // check if this is right
2411                                            kTRUE, // check if this is right
2412                                            AliAODTrack::kPrimary);
2413               vkink->AddDaughter(daughter);
2414               daughter->ConvertAliPIDtoAODPID();
2415             }
2416             else {
2417               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2418               << " track " << idaughter << " has already been used!" << endl;
2419             }
2420           }
2421         }
2422       }
2423     }
2424
2425     // Tracks (primary and orphan)
2426     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2427
2428
2429       if (usedTrack[nTrack]) continue;
2430
2431       AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2432       esdTrack->GetPxPyPz(p);
2433       esdTrack->GetXYZ(pos);
2434       esdTrack->GetCovarianceXYZPxPyPz(covTr);
2435       esdTrack->GetESDpid(pid);
2436
2437       Float_t impactXY, impactZ;
2438
2439       esdTrack->GetImpactParameters(impactXY,impactZ);
2440
2441       if (impactXY<3.) {
2442         // track inside the beam pipe
2443
2444         primary->AddDaughter(aodTrack =
2445             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2446                                          esdTrack->GetLabel(),
2447                                          p,
2448                                          kTRUE,
2449                                          pos,
2450                                          kFALSE,
2451                                          covTr,
2452                                          (Short_t)esdTrack->Charge(),
2453                                          esdTrack->GetITSClusterMap(),
2454                                          pid,
2455                                          primary,
2456                                          kTRUE, // check if this is right
2457                                          kTRUE, // check if this is right
2458                                          AliAODTrack::kPrimary)
2459             );
2460         aodTrack->ConvertAliPIDtoAODPID();
2461       }
2462       else {
2463         // outside the beam pipe: orphan track
2464         // Don't write them anymore!
2465         continue;
2466       }
2467     } // end of loop on tracks
2468
2469     // muon tracks
2470     Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2471     for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2472
2473       AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2474       p[0] = esdMuTrack->Px();
2475       p[1] = esdMuTrack->Py();
2476       p[2] = esdMuTrack->Pz();
2477       pos[0] = primary->GetX();
2478       pos[1] = primary->GetY();
2479       pos[2] = primary->GetZ();
2480
2481       // has to be changed once the muon pid is provided by the ESD
2482       for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2483
2484       primary->AddDaughter(aodTrack =
2485           new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2486                                              0, // no label provided
2487                                              p,
2488                                              kTRUE,
2489                                              pos,
2490                                              kFALSE,
2491                                              NULL, // no covariance matrix provided
2492                                              esdMuTrack->Charge(),
2493                                              0, // ITSClusterMap is set below
2494                                              pid,
2495                                              primary,
2496                                              kFALSE,  // muon tracks are not used to fit the primary vtx
2497                                              kFALSE,  // not used for vertex fit
2498                                              AliAODTrack::kPrimary)
2499           );
2500
2501       aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
2502       Int_t track2Trigger = esdMuTrack->GetMatchTrigger();
2503       aodTrack->SetMatchTrigger(track2Trigger);
2504       if (track2Trigger)
2505         aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
2506       else
2507         aodTrack->SetChi2MatchTrigger(0.);
2508     }
2509
2510     // Access to the AOD container of PMD clusters
2511     TClonesArray &pmdClusters = *(aod->GetPmdClusters());
2512     Int_t jPmdClusters=0;
2513
2514     for (Int_t iPmd = 0; iPmd < nPmdClus; ++iPmd) {
2515       // file pmd clusters, to be revised!
2516       AliESDPmdTrack *pmdTrack = esd->GetPmdTrack(iPmd);
2517       Int_t nLabel = 0;
2518       Int_t *label = 0x0;
2519       Double_t pos[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ() };
2520       Double_t pid[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; // to be revised!
2521       // type not set!
2522       // assoc cluster not set
2523       new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), pos, pid);
2524     }
2525
2526     // Access to the AOD container of clusters
2527     TClonesArray &caloClusters = *(aod->GetCaloClusters());
2528     Int_t jClusters=0;
2529
2530     // Calo Clusters
2531     TArrayS EMCCellNumber(15000);
2532     TArrayD EMCCellAmplitude(15000);
2533     Int_t nEMCCells = 0;
2534     const Float_t fEMCAmpScale = 1./500;
2535
2536     for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
2537
2538       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2539
2540       Int_t id = cluster->GetID();
2541       Int_t nLabel = 0;
2542       Int_t *label = 0x0;
2543       Float_t energy = cluster->E();
2544       cluster->GetPosition(posF);
2545       Char_t ttype=AliAODCluster::kUndef;
2546
2547       if (cluster->GetClusterType() == AliESDCaloCluster::kPHOSCluster) {
2548         ttype=AliAODCluster::kPHOSNeutral;
2549       }
2550       else if (cluster->GetClusterType() == AliESDCaloCluster::kEMCALClusterv1) {
2551         ttype = AliAODCluster::kEMCALClusterv1;
2552       }
2553       else if (cluster->GetClusterType() == AliESDCaloCluster::kEMCALPseudoCluster) {
2554         // Collect raw tower info
2555         for (Int_t iDig = 0; iDig < cluster->GetNumberOfDigits(); iDig++) {
2556           EMCCellNumber[nEMCCells] = cluster->GetDigitIndex()->At(iDig);
2557           EMCCellAmplitude[nEMCCells] = fEMCAmpScale*cluster->GetDigitAmplitude()->At(iDig);
2558           nEMCCells++;
2559         }
2560         // don't write cluster data (it's just a pseudo cluster, holding the tower information)
2561         continue;
2562       }
2563
2564       AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
2565                                                                                         nLabel,
2566                                                                                         label,
2567                                                                                         energy,
2568                                                                                         pos,
2569                                                                                         NULL,
2570                                                                                         ttype);
2571
2572       caloCluster->SetCaloCluster(); // to be refined!
2573
2574     } // end of loop on calo clusters
2575
2576     // fill EMC cell info
2577     AliAODCaloCells &EMCCells = *(aod->GetCaloCells());
2578     EMCCells.CreateContainer(nEMCCells);
2579     EMCCells.SetType(AliAODCaloCells::kEMCAL);
2580     for (Int_t iCell = 0; iCell < nEMCCells; iCell++) {
2581       EMCCells.SetCell(iCell,EMCCellNumber[iCell],EMCCellAmplitude[iCell]);
2582     }
2583     EMCCells.Sort();
2584
2585     // tracklets
2586     AliAODTracklets &SPDTracklets = *(aod->GetTracklets());
2587     const AliMultiplicity *mult = esd->GetMultiplicity();
2588     if (mult) {
2589       if (mult->GetNumberOfTracklets()>0) {
2590         SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
2591
2592         for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
2593           SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n));
2594         }
2595       }
2596     } else {
2597       Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
2598     }
2599
2600     delete [] usedTrack;
2601     delete [] usedV0;
2602     delete [] usedKink;
2603
2604     // fill the tree for this event
2605     aodTree->Fill();
2606   } // end of event loop
2607
2608   aodTree->GetUserInfo()->Add(aod);
2609
2610   // close ESD file
2611   esdFile->Close();
2612
2613   // write the tree to the specified file
2614   aodFile = aodTree->GetCurrentFile();
2615   aodFile->cd();
2616   aodTree->Write();
2617
2618   return;
2619 }
2620
2621 void AliReconstruction::WriteAlignmentData(AliESDEvent* esd)
2622 {
2623   // Write space-points which are then used in the alignment procedures
2624   // For the moment only ITS, TRD and TPC
2625
2626   // Load TOF clusters
2627   if (fTracker[3]){
2628     fLoader[3]->LoadRecPoints("read");
2629     TTree* tree = fLoader[3]->TreeR();
2630     if (!tree) {
2631       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2632       return;
2633     }
2634     fTracker[3]->LoadClusters(tree);
2635   }
2636   Int_t ntracks = esd->GetNumberOfTracks();
2637   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2638     {
2639       AliESDtrack *track = esd->GetTrack(itrack);
2640       Int_t nsp = 0;
2641       Int_t idx[200];
2642       for (Int_t iDet = 3; iDet >= 0; iDet--)
2643         nsp += track->GetNcls(iDet);
2644       if (nsp) {
2645         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2646         track->SetTrackPointArray(sp);
2647         Int_t isptrack = 0;
2648         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2649           AliTracker *tracker = fTracker[iDet];
2650           if (!tracker) continue;
2651           Int_t nspdet = track->GetNcls(iDet);
2652           if (nspdet <= 0) continue;
2653           track->GetClusters(iDet,idx);
2654           AliTrackPoint p;
2655           Int_t isp = 0;
2656           Int_t isp2 = 0;
2657           while (isp < nspdet) {
2658             Bool_t isvalid;
2659             if(IsSelected(fgkDetectorName[iDet],fUseTrackingErrorsForAlignment)) {
2660               isvalid = tracker->GetTrackPointTrackingError(idx[isp2],p,track);
2661             } else {
2662               isvalid = tracker->GetTrackPoint(idx[isp2],p);
2663             }
2664             isp2++;
2665             const Int_t kNTPCmax = 159;
2666             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2667             if (!isvalid) continue;
2668             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2669           }
2670         }
2671       }
2672     }
2673   if (fTracker[3]){
2674     fTracker[3]->UnloadClusters();
2675     fLoader[3]->UnloadRecPoints();
2676   }
2677 }
2678
2679 //_____________________________________________________________________________
2680 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESDEvent* esd)
2681 {
2682   // The method reads the raw-data error log
2683   // accumulated within the rawReader.
2684   // It extracts the raw-data errors related to
2685   // the current event and stores them into
2686   // a TClonesArray inside the esd object.
2687
2688   if (!fRawReader) return;
2689
2690   for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2691
2692     AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2693     if (!log) continue;
2694     if (iEvent != log->GetEventNumber()) continue;
2695
2696     esd->AddRawDataErrorLog(log);
2697   }
2698
2699 }
2700
2701 TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2702   // Dump a file content into a char in TNamed
2703   ifstream in;
2704   in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2705   Int_t kBytes = (Int_t)in.tellg();
2706   printf("Size: %d \n",kBytes);
2707   TNamed *fn = 0;
2708   if(in.good()){
2709     char* memblock = new char [kBytes];
2710     in.seekg (0, ios::beg);
2711     in.read (memblock, kBytes);
2712     in.close();
2713     TString fData(memblock,kBytes);
2714     fn = new TNamed(fName,fData);
2715     printf("fData Size: %d \n",fData.Sizeof());
2716     printf("fName Size: %d \n",fName.Sizeof());
2717     printf("fn    Size: %d \n",fn->Sizeof());
2718     delete[] memblock;
2719   }
2720   else{
2721     AliInfo(Form("Could not Open %s\n",fPath.Data()));
2722   }
2723
2724   return fn;
2725 }
2726
2727 void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2728   // This is not really needed in AliReconstruction at the moment
2729   // but can serve as a template
2730
2731   TList *fList = fTree->GetUserInfo();
2732   TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2733   printf("fn Size: %d \n",fn->Sizeof());
2734
2735   TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2736   const char* cdata = fn->GetTitle();
2737   printf("fTmp Size %d\n",fTmp.Sizeof());
2738
2739   int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2740   printf("calculated size %d\n",size);
2741   ofstream out(fName.Data(),ios::out | ios::binary);
2742   out.write(cdata,size);
2743   out.close();
2744
2745 }
2746
2747
2748
2749 //_____________________________________________________________________________
2750 AliQADataMaker * AliReconstruction::GetQADataMaker(Int_t iDet)
2751 {
2752 // get the quality assurance data maker object and the loader for a detector
2753
2754   if (fQADataMaker[iDet])
2755     return fQADataMaker[iDet];
2756
2757   // load the QA data maker object
2758   TPluginManager* pluginManager = gROOT->GetPluginManager();
2759   TString detName = fgkDetectorName[iDet];
2760   TString qadmName = "Ali" + detName + "QADataMaker";
2761   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT"))
2762     return NULL;
2763
2764   AliQADataMaker * qadm = NULL;
2765   // first check if a plugin is defined for the quality assurance data maker
2766   TPluginHandler* pluginHandler = pluginManager->FindHandler("AliQADataMaker", detName);
2767   // if not, add a plugin for it
2768   if (!pluginHandler) {
2769     AliDebug(1, Form("defining plugin for %s", qadmName.Data()));
2770     TString libs = gSystem->GetLibraries();
2771     if (libs.Contains("lib" + detName + "base.so") ||
2772         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
2773       pluginManager->AddHandler("AliQADataMaker", detName,
2774                                 qadmName, detName + "qadm", qadmName + "()");
2775     } else {
2776       pluginManager->AddHandler("AliQADataMaker", detName,
2777                                 qadmName, detName, qadmName + "()");
2778     }
2779     pluginHandler = pluginManager->FindHandler("AliQADataMaker", detName);
2780   }
2781   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
2782     qadm = (AliQADataMaker *) pluginHandler->ExecPlugin(0);
2783   }
2784   if (qadm) {
2785     AliInfo(Form("Initializing quality assurance data maker for %s", fgkDetectorName[iDet]));
2786     qadm->Init(AliQA::kRECPOINTS, AliCDBManager::Instance()->GetRun(), GetQACycles(fgkDetectorName[iDet]));
2787     qadm->StartOfCycle(AliQA::kRECPOINTS);
2788     qadm->Init(AliQA::kESDS, AliCDBManager::Instance()->GetRun());
2789     qadm->StartOfCycle(AliQA::kESDS, "same") ;
2790     fQADataMaker[iDet] = qadm;
2791   }
2792
2793   return qadm;
2794 }
2795
2796 //_____________________________________________________________________________
2797 Bool_t AliReconstruction::RunQA(const char* detectors, AliESDEvent *& esd)
2798 {
2799   // run the Quality Assurance data producer
2800
2801   AliCodeTimerAuto("")
2802   TString detStr = detectors;
2803   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
2804    if (!IsSelected(fgkDetectorName[iDet], detStr))
2805      continue;
2806    AliQADataMaker * qadm = GetQADataMaker(iDet);
2807    if (!qadm)
2808      continue;
2809    AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2810    AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2811
2812    qadm->Exec(AliQA::kESDS, esd) ;
2813    qadm->Increment() ;
2814
2815    AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2816  }
2817  if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
2818    AliError(Form("the following detectors were not found: %s",
2819                  detStr.Data()));
2820    if (fStopOnError)
2821      return kFALSE;
2822  }
2823
2824  return kTRUE;
2825
2826 }
2827
2828
2829 //_____________________________________________________________________________
2830 void AliReconstruction::CheckQA()
2831 {
2832 // check the QA of SIM for this run and remove the detectors
2833 // with status Fatal
2834
2835         TString newDetList ;
2836         for (Int_t iDet = 0; iDet < AliQA::kNDET; iDet++) {
2837                 TString detName(AliQA::GetDetName(iDet)) ;
2838                 if ( fRunLocalReconstruction.Contains(AliQA::GetDetName(iDet)) ||
2839                         fRunLocalReconstruction.Contains("ALL") )  {
2840                         AliQA * qa = AliQA::Instance(AliQA::DETECTORINDEX(iDet)) ;
2841                         if ( qa->IsSet(AliQA::DETECTORINDEX(iDet), AliQA::kSIM, AliQA::kFATAL)) {
2842                                 AliInfo(Form("QA status for %s in Hits and/or SDIGITS  and/or Digits was Fatal; No reconstruction performed", detName.Data())) ;
2843                         } else if ( qa->IsSet(AliQA::DETECTORINDEX(iDet), AliQA::kSIM, AliQA::kERROR)) {
2844                                 AliError(Form("QA status for %s in Hits and/or SDIGITS  and/or Digits was ERROR", detName.Data())) ;
2845                                 newDetList += detName ;
2846                                 newDetList += " " ;
2847                         } else if ( qa->IsSet(AliQA::DETECTORINDEX(iDet), AliQA::kSIM, AliQA::kWARNING) ) {
2848                                 AliWarning(Form("QA status for %s in Hits and/or SDIGITS  and/or Digits was WARNING", detName.Data())) ;
2849                                 newDetList += detName ;
2850                                 newDetList += " " ;
2851                         } else if ( qa->IsSet(AliQA::DETECTORINDEX(iDet), AliQA::kSIM, AliQA::kINFO) ) {
2852                                 AliInfo(Form("QA status for %s in Hits and/or SDIGITS  and/or Digits was INFO", detName.Data())) ;
2853                                 newDetList += detName ;
2854                                 newDetList += " " ;
2855                         }
2856                 }
2857         }
2858         fRunLocalReconstruction = newDetList ;
2859 }
2860
2861 //_____________________________________________________________________________
2862 Int_t AliReconstruction::GetDetIndex(const char* detector)
2863 {
2864   // return the detector index corresponding to detector
2865   Int_t index = -1 ;
2866   for (index = 0; index < fgkNDetectors ; index++) {
2867     if ( strcmp(detector, fgkDetectorName[index]) == 0 )
2868         break ;
2869   }
2870   return index ;
2871 }