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