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