STEER/AliReconstruction.cxx

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