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