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 // The name of the galice file can be changed from the default               //
  51 // "galice.root" by passing it as argument to the AliReconstruction          //
  52 // constructor or by                                                         //
  53 //                                                                           //
  54 //   rec.SetGAliceFile("...");                                               //
  55 //                                                                           //
  56 // The local reconstruction can be switched on or off for individual         //
  57 // detectors by                                                              //
  58 //                                                                           //
  59 //   rec.SetRunLocalReconstruction("...");                                   //
  60 //                                                                           //
  61 // The argument is a (case sensitive) string with the names of the           //
  62 // detectors separated by a space. The special string "ALL" selects all      //
  63 // available detectors. This is the default.                                 //
  64 //                                                                           //
  65 // The reconstruction of the primary vertex position can be switched off by  //
  66 //                                                                           //
  67 //   rec.SetRunVertexFinder(kFALSE);                                         //
  68 //                                                                           //
  69 // The tracking and the creation of ESD tracks can be switched on for        //
  70 // selected detectors by                                                     //
  71 //                                                                           //
  72 //   rec.SetRunTracking("...");                                              //
  73 //                                                                           //
  74 // Uniform/nonuniform field tracking switches (default: uniform field)       //
  75 //                                                                           //
  76 //   rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
  77 //                                                                           //
  78 // The filling of additional ESD information can be steered by               //
  79 //                                                                           //
  80 //   rec.SetFillESD("...");                                                  //
  81 //                                                                           //
  82 // Again, for both methods the string specifies the list of detectors.       //
  83 // The default is "ALL".                                                     //
  84 //                                                                           //
  85 // The call of the shortcut method                                           //
  86 //                                                                           //
  87 //   rec.SetRunReconstruction("...");                                        //
  88 //                                                                           //
  89 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
  90 // SetFillESD with the same detector selecting string as argument.           //
  91 //                                                                           //
  92 // The reconstruction requires digits or raw data as input. For the creation //
  93 // of digits and raw data have a look at the class AliSimulation.            //
  94 //                                                                           //
  95 // For debug purposes the method SetCheckPointLevel can be used. If the      //
  96 // argument is greater than 0, files with ESD events will be written after   //
  97 // selected steps of the reconstruction for each event:                      //
  98 //   level 1: after tracking and after filling of ESD (final)                //
  99 //   level 2: in addition after each tracking step                           //
 100 //   level 3: in addition after the filling of ESD for each detector         //
 101 // If a final check point file exists for an event, this event will be       //
 102 // skipped in the reconstruction. The tracking and the filling of ESD for    //
 103 // a detector will be skipped as well, if the corresponding check point      //
 104 // file exists. The ESD event will then be loaded from the file instead.     //
 105 //                                                                           //
 106 ///////////////////////////////////////////////////////////////////////////////
 107
 108 #include <TArrayF.h>
 109 #include <TFile.h>
 110 #include <TSystem.h>
 111 #include <TROOT.h>
 112 #include <TPluginManager.h>
 113 #include <TStopwatch.h>
 114 #include <TGeoManager.h>
 115 #include <TLorentzVector.h>
 116
 117 #include "AliReconstruction.h"
 118 #include "AliReconstructor.h"
 119 #include "AliLog.h"
 120 #include "AliRunLoader.h"
 121 #include "AliRun.h"
 122 #include "AliRawReaderFile.h"
 123 #include "AliRawReaderDate.h"
 124 #include "AliRawReaderRoot.h"
 125 #include "AliRawEventHeaderBase.h"
 126 #include "AliESD.h"
 127 #include "AliESDfriend.h"
 128 #include "AliESDVertex.h"
 129 #include "AliMultiplicity.h"
 130 #include "AliTracker.h"
 131 #include "AliVertexer.h"
 132 #include "AliVertexerTracks.h"
 133 #include "AliHeader.h"
 134 #include "AliGenEventHeader.h"
 135 #include "AliPID.h"
 136 #include "AliESDpid.h"
 137 #include "AliESDtrack.h"
 138
 139 #include "AliRunTag.h"
 140 #include "AliDetectorTag.h"
 141 #include "AliEventTag.h"
 142
 143 #include "AliTrackPointArray.h"
 144 #include "AliCDBManager.h"
 145 #include "AliCDBEntry.h"
 146 #include "AliAlignObj.h"
 147
 148 #include "AliCentralTrigger.h"
 149 #include "AliCTPRawStream.h"
 150
 151 ClassImp(AliReconstruction)
 152
 153
 154 //_____________________________________________________________________________
 155 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
 156
 157 //_____________________________________________________________________________
 158 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
 159                                      const char* name, const char* title) :
 160   TNamed(name, title),
 161
 162   fUniformField(kTRUE),
 163   fRunVertexFinder(kTRUE),
 164   fRunHLTTracking(kFALSE),
 165   fStopOnError(kFALSE),
 166   fWriteAlignmentData(kFALSE),
 167   fWriteESDfriend(kFALSE),
 168   fFillTriggerESD(kTRUE),
 169
 170   fRunLocalReconstruction("ALL"),
 171   fRunTracking("ALL"),
 172   fFillESD("ALL"),
 173   fGAliceFileName(gAliceFilename),
 174   fInput(""),
 175   fEquipIdMap(""),
 176   fFirstEvent(0),
 177   fLastEvent(-1),
 178   fCheckPointLevel(0),
 179   fOptions(),
 180   fLoadAlignFromCDB(kTRUE),
 181   fLoadAlignData("ALL"),
 182
 183   fRunLoader(NULL),
 184   fRawReader(NULL),
 185
 186   fVertexer(NULL),
 187   fDiamondProfile(NULL),
 188
 189   fAlignObjArray(NULL),
 190   fCDBUri(cdbUri),
 191   fSpecCDBUri()
 192 {
 193 // create reconstruction object with default parameters
 194
 195   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 196     fReconstructor[iDet] = NULL;
 197     fLoader[iDet] = NULL;
 198     fTracker[iDet] = NULL;
 199   }
 200   AliPID pid;
 201 }
 202
 203 //_____________________________________________________________________________
 204 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
 205   TNamed(rec),
 206
 207   fUniformField(rec.fUniformField),
 208   fRunVertexFinder(rec.fRunVertexFinder),
 209   fRunHLTTracking(rec.fRunHLTTracking),
 210   fStopOnError(rec.fStopOnError),
 211   fWriteAlignmentData(rec.fWriteAlignmentData),
 212   fWriteESDfriend(rec.fWriteESDfriend),
 213   fFillTriggerESD(rec.fFillTriggerESD),
 214
 215   fRunLocalReconstruction(rec.fRunLocalReconstruction),
 216   fRunTracking(rec.fRunTracking),
 217   fFillESD(rec.fFillESD),
 218   fGAliceFileName(rec.fGAliceFileName),
 219   fInput(rec.fInput),
 220   fEquipIdMap(rec.fEquipIdMap),
 221   fFirstEvent(rec.fFirstEvent),
 222   fLastEvent(rec.fLastEvent),
 223   fCheckPointLevel(0),
 224   fOptions(),
 225   fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
 226   fLoadAlignData(rec.fLoadAlignData),
 227
 228   fRunLoader(NULL),
 229   fRawReader(NULL),
 230
 231   fVertexer(NULL),
 232   fDiamondProfile(NULL),
 233
 234   fAlignObjArray(rec.fAlignObjArray),
 235   fCDBUri(rec.fCDBUri),
 236   fSpecCDBUri()
 237 {
 238 // copy constructor
 239
 240   for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
 241     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
 242   }
 243   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 244     fReconstructor[iDet] = NULL;
 245     fLoader[iDet] = NULL;
 246     fTracker[iDet] = NULL;
 247   }
 248   for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
 249     if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
 250   }
 251 }
 252
 253 //_____________________________________________________________________________
 254 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 255 {
 256 // assignment operator
 257
 258   this->~AliReconstruction();
 259   new(this) AliReconstruction(rec);
 260   return *this;
 261 }
 262
 263 //_____________________________________________________________________________
 264 AliReconstruction::~AliReconstruction()
 265 {
 266 // clean up
 267
 268   CleanUp();
 269   fOptions.Delete();
 270   fSpecCDBUri.Delete();
 271 }
 272
 273 //_____________________________________________________________________________
 274 void AliReconstruction::InitCDBStorage()
 275 {
 276 // activate a default CDB storage
 277 // First check if we have any CDB storage set, because it is used
 278 // to retrieve the calibration and alignment constants
 279
 280   AliCDBManager* man = AliCDBManager::Instance();
 281   if (man->IsDefaultStorageSet())
 282   {
 283     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 284     AliWarning("Default CDB storage has been already set !");
 285     AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
 286     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 287     fCDBUri = "";
 288   }
 289   else {
 290     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 291     AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
 292     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 293     man->SetDefaultStorage(fCDBUri);
 294   }
 295
 296   // Now activate the detector specific CDB storage locations
 297   for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
 298     TObject* obj = fSpecCDBUri[i];
 299     if (!obj) continue;
 300     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 301     AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
 302     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 303     man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
 304   }
 305   man->Print();
 306 }
 307
 308 //_____________________________________________________________________________
 309 void AliReconstruction::SetDefaultStorage(const char* uri) {
 310 // Store the desired default CDB storage location
 311 // Activate it later within the Run() method
 312
 313   fCDBUri = uri;
 314
 315 }
 316
 317 //_____________________________________________________________________________
 318 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
 319 // Store a detector-specific CDB storage location
 320 // Activate it later within the Run() method
 321
 322   AliCDBPath aPath(calibType);
 323   if(!aPath.IsValid()){
 324         // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
 325         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 326                 if(!strcmp(calibType, fgkDetectorName[iDet])) {
 327                         aPath.SetPath(Form("%s/*", calibType));
 328                         AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
 329                         break;
 330                 }
 331         }
 332         if(!aPath.IsValid()){
 333                 AliError(Form("Not a valid path or detector: %s", calibType));
 334                 return;
 335         }
 336   }
 337
 338   // check that calibType refers to a "valid" detector name
 339   Bool_t isDetector = kFALSE;
 340   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 341     TString detName = fgkDetectorName[iDet];
 342     if(aPath.GetLevel0() == detName) {
 343         isDetector = kTRUE;
 344         break;
 345     }
 346   }
 347
 348   if(!isDetector) {
 349         AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
 350         return;
 351   }
 352
 353   TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
 354   if (obj) fSpecCDBUri.Remove(obj);
 355   fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
 356
 357 }
 358
 359 //_____________________________________________________________________________
 360 Bool_t AliReconstruction::SetRunNumber()
 361 {
 362   // The method is called in Run() in order
 363   // to set a correct run number.
 364   // In case of raw data reconstruction the
 365   // run number is taken from the raw data header
 366
 367   if(AliCDBManager::Instance()->GetRun() < 0) {
 368     if (!fRunLoader) {
 369       AliError("No run loader is found !");
 370       return kFALSE;
 371     }
 372     // read run number from gAlice
 373     if(fRunLoader->GetAliRun())
 374       AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
 375     else {
 376       if(fRawReader) {
 377         if(fRawReader->NextEvent()) {
 378           AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
 379           fRawReader->RewindEvents();
 380         }
 381         else {
 382           AliError("No raw-data events found !");
 383           return kFALSE;
 384         }
 385       }
 386       else {
 387         AliError("Neither gAlice nor RawReader objects are found !");
 388         return kFALSE;
 389       }
 390     }
 391     AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
 392   }
 393   return kTRUE;
 394 }
 395
 396 //_____________________________________________________________________________
 397 Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
 398 {
 399   // Read collection of alignment objects (AliAlignObj derived) saved
 400   // in the TClonesArray ClArrayName and apply them to the geometry
 401   // manager singleton.
 402   //
 403   alObjArray->Sort();
 404   Int_t nvols = alObjArray->GetEntriesFast();
 405
 406   Bool_t flag = kTRUE;
 407
 408   for(Int_t j=0; j<nvols; j++)
 409     {
 410       AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
 411       if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
 412     }
 413
 414   if (AliDebugLevelClass() >= 1) {
 415     gGeoManager->GetTopNode()->CheckOverlaps(20);
 416     TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
 417     if(ovexlist->GetEntriesFast()){
 418       AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
 419    }
 420   }
 421
 422   return flag;
 423
 424 }
 425
 426 //_____________________________________________________________________________
 427 Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
 428 {
 429   // Fills array of single detector's alignable objects from CDB
 430
 431   AliDebug(2, Form("Loading alignment data for detector: %s",detName));
 432
 433   AliCDBEntry *entry;
 434
 435   AliCDBPath path(detName,"Align","Data");
 436
 437   entry=AliCDBManager::Instance()->Get(path.GetPath());
 438   if(!entry){
 439         AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
 440         return kFALSE;
 441   }
 442   entry->SetOwner(1);
 443   TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
 444   alignArray->SetOwner(0);
 445   AliDebug(2,Form("Found %d alignment objects for %s",
 446                         alignArray->GetEntries(),detName));
 447
 448   AliAlignObj *alignObj=0;
 449   TIter iter(alignArray);
 450
 451   // loop over align objects in detector
 452   while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
 453         fAlignObjArray->Add(alignObj);
 454   }
 455   // delete entry --- Don't delete, it is cached!
 456
 457   AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
 458   return kTRUE;
 459
 460 }
 461
 462 //_____________________________________________________________________________
 463 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
 464 {
 465   // Read the alignment objects from CDB.
 466   // Each detector is supposed to have the
 467   // alignment objects in DET/Align/Data CDB path.
 468   // All the detector objects are then collected,
 469   // sorted by geometry level (starting from ALIC) and
 470   // then applied to the TGeo geometry.
 471   // Finally an overlaps check is performed.
 472
 473   // Load alignment data from CDB and fill fAlignObjArray
 474   if(fLoadAlignFromCDB){
 475         if(!fAlignObjArray) fAlignObjArray = new TObjArray();
 476
 477         //fAlignObjArray->RemoveAll();
 478         fAlignObjArray->Clear();
 479         fAlignObjArray->SetOwner(0);
 480
 481         TString detStr = detectors;
 482         TString dataNotLoaded="";
 483         TString dataLoaded="";
 484
 485         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 486           if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 487           if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
 488             dataNotLoaded += fgkDetectorName[iDet];
 489             dataNotLoaded += " ";
 490           } else {
 491             dataLoaded += fgkDetectorName[iDet];
 492             dataLoaded += " ";
 493           }
 494         } // end loop over detectors
 495
 496         if ((detStr.CompareTo("ALL") == 0)) detStr = "";
 497         dataNotLoaded += detStr;
 498         if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
 499                           dataLoaded.Data()));
 500         if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
 501                           dataNotLoaded.Data()));
 502   } // fLoadAlignFromCDB flag
 503
 504   // Check if the array with alignment objects was
 505   // provided by the user. If yes, apply the objects
 506   // to the present TGeo geometry
 507   if (fAlignObjArray) {
 508     if (gGeoManager && gGeoManager->IsClosed()) {
 509       if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
 510         AliError("The misalignment of one or more volumes failed!"
 511                  "Compare the list of simulated detectors and the list of detector alignment data!");
 512         return kFALSE;
 513       }
 514     }
 515     else {
 516       AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
 517       return kFALSE;
 518     }
 519   }
 520
 521   return kTRUE;
 522 }
 523
 524 //_____________________________________________________________________________
 525 void AliReconstruction::SetGAliceFile(const char* fileName)
 526 {
 527 // set the name of the galice file
 528
 529   fGAliceFileName = fileName;
 530 }
 531
 532 //_____________________________________________________________________________
 533 void AliReconstruction::SetOption(const char* detector, const char* option)
 534 {
 535 // set options for the reconstruction of a detector
 536
 537   TObject* obj = fOptions.FindObject(detector);
 538   if (obj) fOptions.Remove(obj);
 539   fOptions.Add(new TNamed(detector, option));
 540 }
 541
 542
 543 //_____________________________________________________________________________
 544 Bool_t AliReconstruction::Run(const char* input,
 545                               Int_t firstEvent, Int_t lastEvent)
 546 {
 547 // run the reconstruction
 548
 549   // set the input
 550   if (!input) input = fInput.Data();
 551   TString fileName(input);
 552   if (fileName.EndsWith("/")) {
 553     fRawReader = new AliRawReaderFile(fileName);
 554   } else if (fileName.EndsWith(".root")) {
 555     fRawReader = new AliRawReaderRoot(fileName);
 556   } else if (!fileName.IsNull()) {
 557     fRawReader = new AliRawReaderDate(fileName);
 558     fRawReader->SelectEvents(7);
 559   }
 560   if (!fEquipIdMap.IsNull() && fRawReader)
 561     fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
 562
 563
 564   // get the run loader
 565   if (!InitRunLoader()) return kFALSE;
 566
 567   // Initialize the CDB storage
 568   InitCDBStorage();
 569
 570   // Set run number in CDBManager (if it is not already set by the user)
 571   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
 572
 573   // Import ideal TGeo geometry and apply misalignment
 574   if (!gGeoManager) {
 575     TString geom(gSystem->DirName(fGAliceFileName));
 576     geom += "/geometry.root";
 577     TGeoManager::Import(geom.Data());
 578     if (!gGeoManager) if (fStopOnError) return kFALSE;
 579   }
 580   if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
 581
 582   // local reconstruction
 583   if (!fRunLocalReconstruction.IsNull()) {
 584     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
 585       if (fStopOnError) {CleanUp(); return kFALSE;}
 586     }
 587   }
 588 //  if (!fRunVertexFinder && fRunTracking.IsNull() &&
 589 //      fFillESD.IsNull()) return kTRUE;
 590
 591   // get vertexer
 592   if (fRunVertexFinder && !CreateVertexer()) {
 593     if (fStopOnError) {
 594       CleanUp();
 595       return kFALSE;
 596     }
 597   }
 598
 599   // get trackers
 600   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
 601     if (fStopOnError) {
 602       CleanUp();
 603       return kFALSE;
 604     }
 605   }
 606
 607
 608   TStopwatch stopwatch;
 609   stopwatch.Start();
 610
 611   // get the possibly already existing ESD file and tree
 612   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
 613   TFile* fileOld = NULL;
 614   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
 615   if (!gSystem->AccessPathName("AliESDs.root")){
 616     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
 617     fileOld = TFile::Open("AliESDs.old.root");
 618     if (fileOld && fileOld->IsOpen()) {
 619       treeOld = (TTree*) fileOld->Get("esdTree");
 620       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
 621       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
 622       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
 623     }
 624   }
 625
 626   // create the ESD output file and tree
 627   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 628   if (!file->IsOpen()) {
 629     AliError("opening AliESDs.root failed");
 630     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 631   }
 632   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 633   tree->Branch("ESD", "AliESD", &esd);
 634   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
 635   hlttree->Branch("ESD", "AliESD", &hltesd);
 636   delete esd; delete hltesd;
 637   esd = NULL; hltesd = NULL;
 638
 639   // create the branch with ESD additions
 640   AliESDfriend *esdf=0;
 641   if (fWriteESDfriend) {
 642      TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
 643      br->SetFile("AliESDfriends.root");
 644   }
 645
 646   AliVertexerTracks tVertexer;
 647   if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
 648
 649   // loop over events
 650   if (fRawReader) fRawReader->RewindEvents();
 651
 652   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 653     if (fRawReader) fRawReader->NextEvent();
 654     if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
 655       // copy old ESD to the new one
 656       if (treeOld) {
 657         treeOld->SetBranchAddress("ESD", &esd);
 658         treeOld->GetEntry(iEvent);
 659       }
 660       tree->Fill();
 661       if (hlttreeOld) {
 662         hlttreeOld->SetBranchAddress("ESD", &hltesd);
 663         hlttreeOld->GetEntry(iEvent);
 664       }
 665       hlttree->Fill();
 666       continue;
 667     }
 668
 669     AliInfo(Form("processing event %d", iEvent));
 670     fRunLoader->GetEvent(iEvent);
 671
 672     char fileName[256];
 673     sprintf(fileName, "ESD_%d.%d_final.root",
 674             fRunLoader->GetHeader()->GetRun(),
 675             fRunLoader->GetHeader()->GetEventNrInRun());
 676     if (!gSystem->AccessPathName(fileName)) continue;
 677
 678     // local reconstruction
 679     if (!fRunLocalReconstruction.IsNull()) {
 680       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
 681         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 682       }
 683     }
 684
 685     esd = new AliESD; hltesd = new AliESD;
 686     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 687     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 688     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 689     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 690
 691     // Set magnetic field from the tracker
 692     esd->SetMagneticField(AliTracker::GetBz());
 693     hltesd->SetMagneticField(AliTracker::GetBz());
 694
 695     // vertex finder
 696     if (fRunVertexFinder) {
 697       if (!ReadESD(esd, "vertex")) {
 698         if (!RunVertexFinder(esd)) {
 699           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 700         }
 701         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 702       }
 703     }
 704
 705     // HLT tracking
 706     if (!fRunTracking.IsNull()) {
 707       if (fRunHLTTracking) {
 708         hltesd->SetVertex(esd->GetVertex());
 709         if (!RunHLTTracking(hltesd)) {
 710           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 711         }
 712       }
 713     }
 714
 715     // barrel tracking
 716     if (!fRunTracking.IsNull()) {
 717       if (!ReadESD(esd, "tracking")) {
 718         if (!RunTracking(esd)) {
 719           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 720         }
 721         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 722       }
 723     }
 724
 725     // fill ESD
 726     if (!fFillESD.IsNull()) {
 727       if (!FillESD(esd, fFillESD)) {
 728         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 729       }
 730     }
 731     // fill Event header information from the RawEventHeader
 732     if (fRawReader){FillRawEventHeaderESD(esd);}
 733
 734     // combined PID
 735     AliESDpid::MakePID(esd);
 736     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 737
 738     if (fFillTriggerESD) {
 739       if (!ReadESD(esd, "trigger")) {
 740         if (!FillTriggerESD(esd)) {
 741           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 742         }
 743         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
 744       }
 745     }
 746
 747     esd->SetPrimaryVertex(tVertexer.FindPrimaryVertex(esd));
 748
 749     // write ESD
 750     if (fWriteESDfriend) {
 751        esdf=new AliESDfriend();
 752        esd->GetESDfriend(esdf);
 753     }
 754     tree->Fill();
 755
 756     // write HLT ESD
 757     hlttree->Fill();
 758
 759     if (fCheckPointLevel > 0)  WriteESD(esd, "final");
 760
 761     delete esd; delete esdf; delete hltesd;
 762     esd = NULL; esdf=NULL; hltesd = NULL;
 763   }
 764
 765   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
 766                stopwatch.RealTime(),stopwatch.CpuTime()));
 767
 768   file->cd();
 769   if (fWriteESDfriend)
 770     tree->SetBranchStatus("ESDfriend*",0);
 771   tree->Write();
 772   hlttree->Write();
 773
 774   // Create tags for the events in the ESD tree (the ESD tree is always present)
 775   // In case of empty events the tags will contain dummy values
 776   CreateTag(file);
 777   CleanUp(file, fileOld);
 778
 779   return kTRUE;
 780 }
 781
 782
 783 //_____________________________________________________________________________
 784 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 785 {
 786 // run the local reconstruction
 787
 788   TStopwatch stopwatch;
 789   stopwatch.Start();
 790
 791   TString detStr = detectors;
 792   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 793     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 794     AliReconstructor* reconstructor = GetReconstructor(iDet);
 795     if (!reconstructor) continue;
 796     if (reconstructor->HasLocalReconstruction()) continue;
 797
 798     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 799     TStopwatch stopwatchDet;
 800     stopwatchDet.Start();
 801     if (fRawReader) {
 802       fRawReader->RewindEvents();
 803       reconstructor->Reconstruct(fRunLoader, fRawReader);
 804     } else {
 805       reconstructor->Reconstruct(fRunLoader);
 806     }
 807     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 808                  fgkDetectorName[iDet],
 809                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 810   }
 811
 812   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 813     AliError(Form("the following detectors were not found: %s",
 814                   detStr.Data()));
 815     if (fStopOnError) return kFALSE;
 816   }
 817
 818   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 819                stopwatch.RealTime(),stopwatch.CpuTime()));
 820
 821   return kTRUE;
 822 }
 823
 824 //_____________________________________________________________________________
 825 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 826 {
 827 // run the local reconstruction
 828
 829   TStopwatch stopwatch;
 830   stopwatch.Start();
 831
 832   TString detStr = detectors;
 833   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 834     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 835     AliReconstructor* reconstructor = GetReconstructor(iDet);
 836     if (!reconstructor) continue;
 837     AliLoader* loader = fLoader[iDet];
 838
 839     // conversion of digits
 840     if (fRawReader && reconstructor->HasDigitConversion()) {
 841       AliInfo(Form("converting raw data digits into root objects for %s",
 842                    fgkDetectorName[iDet]));
 843       TStopwatch stopwatchDet;
 844       stopwatchDet.Start();
 845       loader->LoadDigits("update");
 846       loader->CleanDigits();
 847       loader->MakeDigitsContainer();
 848       TTree* digitsTree = loader->TreeD();
 849       reconstructor->ConvertDigits(fRawReader, digitsTree);
 850       loader->WriteDigits("OVERWRITE");
 851       loader->UnloadDigits();
 852       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 853                    fgkDetectorName[iDet],
 854                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 855     }
 856
 857     // local reconstruction
 858     if (!reconstructor->HasLocalReconstruction()) continue;
 859     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 860     TStopwatch stopwatchDet;
 861     stopwatchDet.Start();
 862     loader->LoadRecPoints("update");
 863     loader->CleanRecPoints();
 864     loader->MakeRecPointsContainer();
 865     TTree* clustersTree = loader->TreeR();
 866     if (fRawReader && !reconstructor->HasDigitConversion()) {
 867       reconstructor->Reconstruct(fRawReader, clustersTree);
 868     } else {
 869       loader->LoadDigits("read");
 870       TTree* digitsTree = loader->TreeD();
 871       if (!digitsTree) {
 872         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 873         if (fStopOnError) return kFALSE;
 874       } else {
 875         reconstructor->Reconstruct(digitsTree, clustersTree);
 876       }
 877       loader->UnloadDigits();
 878     }
 879     loader->WriteRecPoints("OVERWRITE");
 880     loader->UnloadRecPoints();
 881     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
 882                     fgkDetectorName[iDet],
 883                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 884   }
 885
 886   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 887     AliError(Form("the following detectors were not found: %s",
 888                   detStr.Data()));
 889     if (fStopOnError) return kFALSE;
 890   }
 891
 892   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 893                stopwatch.RealTime(),stopwatch.CpuTime()));
 894
 895   return kTRUE;
 896 }
 897
 898 //_____________________________________________________________________________
 899 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 900 {
 901 // run the barrel tracking
 902
 903   TStopwatch stopwatch;
 904   stopwatch.Start();
 905
 906   AliESDVertex* vertex = NULL;
 907   Double_t vtxPos[3] = {0, 0, 0};
 908   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 909   TArrayF mcVertex(3);
 910   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 911     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 912     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 913   }
 914
 915   if (fVertexer) {
 916     AliInfo("running the ITS vertex finder");
 917     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 918     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 919     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 920     if(!vertex){
 921       AliWarning("Vertex not found");
 922       vertex = new AliESDVertex();
 923       vertex->SetName("default");
 924     }
 925     else {
 926       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 927       vertex->SetName("reconstructed");
 928     }
 929
 930   } else {
 931     AliInfo("getting the primary vertex from MC");
 932     vertex = new AliESDVertex(vtxPos, vtxErr);
 933   }
 934
 935   if (vertex) {
 936     vertex->GetXYZ(vtxPos);
 937     vertex->GetSigmaXYZ(vtxErr);
 938   } else {
 939     AliWarning("no vertex reconstructed");
 940     vertex = new AliESDVertex(vtxPos, vtxErr);
 941   }
 942   esd->SetVertex(vertex);
 943   // if SPD multiplicity has been determined, it is stored in the ESD
 944   AliMultiplicity *mult= fVertexer->GetMultiplicity();
 945   if(mult)esd->SetMultiplicity(mult);
 946
 947   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 948     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
 949   }
 950   delete vertex;
 951
 952   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 953                stopwatch.RealTime(),stopwatch.CpuTime()));
 954
 955   return kTRUE;
 956 }
 957
 958 //_____________________________________________________________________________
 959 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
 960 {
 961 // run the HLT barrel tracking
 962
 963   TStopwatch stopwatch;
 964   stopwatch.Start();
 965
 966   if (!fRunLoader) {
 967     AliError("Missing runLoader!");
 968     return kFALSE;
 969   }
 970
 971   AliInfo("running HLT tracking");
 972
 973   // Get a pointer to the HLT reconstructor
 974   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
 975   if (!reconstructor) return kFALSE;
 976
 977   // TPC + ITS
 978   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 979     TString detName = fgkDetectorName[iDet];
 980     AliDebug(1, Form("%s HLT tracking", detName.Data()));
 981     reconstructor->SetOption(detName.Data());
 982     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
 983     if (!tracker) {
 984       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
 985       if (fStopOnError) return kFALSE;
 986       continue;
 987     }
 988     Double_t vtxPos[3];
 989     Double_t vtxErr[3]={0.005,0.005,0.010};
 990     const AliESDVertex *vertex = esd->GetVertex();
 991     vertex->GetXYZ(vtxPos);
 992     tracker->SetVertex(vtxPos,vtxErr);
 993     if(iDet != 1) {
 994       fLoader[iDet]->LoadRecPoints("read");
 995       TTree* tree = fLoader[iDet]->TreeR();
 996       if (!tree) {
 997         AliError(Form("Can't get the %s cluster tree", detName.Data()));
 998         return kFALSE;
 999       }
1000       tracker->LoadClusters(tree);
1001     }
1002     if (tracker->Clusters2Tracks(esd) != 0) {
1003       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1004       return kFALSE;
1005     }
1006     if(iDet != 1) {
1007       tracker->UnloadClusters();
1008     }
1009     delete tracker;
1010   }
1011
1012   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1013                stopwatch.RealTime(),stopwatch.CpuTime()));
1014
1015   return kTRUE;
1016 }
1017
1018 //_____________________________________________________________________________
1019 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1020 {
1021 // run the barrel tracking
1022
1023   TStopwatch stopwatch;
1024   stopwatch.Start();
1025
1026   AliInfo("running tracking");
1027
1028   // pass 1: TPC + ITS inwards
1029   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1030     if (!fTracker[iDet]) continue;
1031     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1032
1033     // load clusters
1034     fLoader[iDet]->LoadRecPoints("read");
1035     TTree* tree = fLoader[iDet]->TreeR();
1036     if (!tree) {
1037       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1038       return kFALSE;
1039     }
1040     fTracker[iDet]->LoadClusters(tree);
1041
1042     // run tracking
1043     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1044       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1045       return kFALSE;
1046     }
1047     if (fCheckPointLevel > 1) {
1048       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1049     }
1050     // preliminary PID in TPC needed by the ITS tracker
1051     if (iDet == 1) {
1052       GetReconstructor(1)->FillESD(fRunLoader, esd);
1053       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1054       AliESDpid::MakePID(esd);
1055     }
1056   }
1057
1058   // pass 2: ALL backwards
1059   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1060     if (!fTracker[iDet]) continue;
1061     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1062
1063     // load clusters
1064     if (iDet > 1) {     // all except ITS, TPC
1065       TTree* tree = NULL;
1066       fLoader[iDet]->LoadRecPoints("read");
1067       tree = fLoader[iDet]->TreeR();
1068       if (!tree) {
1069         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1070         return kFALSE;
1071       }
1072       fTracker[iDet]->LoadClusters(tree);
1073     }
1074
1075     // run tracking
1076     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1077       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1078       return kFALSE;
1079     }
1080     if (fCheckPointLevel > 1) {
1081       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1082     }
1083
1084     // unload clusters
1085     if (iDet > 2) {     // all except ITS, TPC, TRD
1086       fTracker[iDet]->UnloadClusters();
1087       fLoader[iDet]->UnloadRecPoints();
1088     }
1089     // updated PID in TPC needed by the ITS tracker -MI
1090     if (iDet == 1) {
1091       GetReconstructor(1)->FillESD(fRunLoader, esd);
1092       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1093       AliESDpid::MakePID(esd);
1094     }
1095   }
1096
1097   // write space-points to the ESD in case alignment data output
1098   // is switched on
1099   if (fWriteAlignmentData)
1100     WriteAlignmentData(esd);
1101
1102   // pass 3: TRD + TPC + ITS refit inwards
1103   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1104     if (!fTracker[iDet]) continue;
1105     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1106
1107     // run tracking
1108     if (fTracker[iDet]->RefitInward(esd) != 0) {
1109       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1110       return kFALSE;
1111     }
1112     if (fCheckPointLevel > 1) {
1113       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1114     }
1115
1116     // unload clusters
1117     fTracker[iDet]->UnloadClusters();
1118     fLoader[iDet]->UnloadRecPoints();
1119   }
1120   //
1121   // Propagate track to the vertex - if not done by ITS
1122   //
1123   Int_t ntracks = esd->GetNumberOfTracks();
1124   for (Int_t itrack=0; itrack<ntracks; itrack++){
1125     const Double_t kRadius  = 3;   // beam pipe radius
1126     const Double_t kMaxStep = 5;   // max step
1127     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1128     Double_t       fieldZ   = AliTracker::GetBz();  //
1129     AliESDtrack * track = esd->GetTrack(itrack);
1130     if (!track) continue;
1131     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1132     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1133     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1134   }
1135
1136   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1137                stopwatch.RealTime(),stopwatch.CpuTime()));
1138
1139   return kTRUE;
1140 }
1141
1142 //_____________________________________________________________________________
1143 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1144 {
1145 // fill the event summary data
1146
1147   TStopwatch stopwatch;
1148   stopwatch.Start();
1149   AliInfo("filling ESD");
1150
1151   TString detStr = detectors;
1152   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1153     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1154     AliReconstructor* reconstructor = GetReconstructor(iDet);
1155     if (!reconstructor) continue;
1156
1157     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1158       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1159       TTree* clustersTree = NULL;
1160       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1161         fLoader[iDet]->LoadRecPoints("read");
1162         clustersTree = fLoader[iDet]->TreeR();
1163         if (!clustersTree) {
1164           AliError(Form("Can't get the %s clusters tree",
1165                         fgkDetectorName[iDet]));
1166           if (fStopOnError) return kFALSE;
1167         }
1168       }
1169       if (fRawReader && !reconstructor->HasDigitConversion()) {
1170         reconstructor->FillESD(fRawReader, clustersTree, esd);
1171       } else {
1172         TTree* digitsTree = NULL;
1173         if (fLoader[iDet]) {
1174           fLoader[iDet]->LoadDigits("read");
1175           digitsTree = fLoader[iDet]->TreeD();
1176           if (!digitsTree) {
1177             AliError(Form("Can't get the %s digits tree",
1178                           fgkDetectorName[iDet]));
1179             if (fStopOnError) return kFALSE;
1180           }
1181         }
1182         reconstructor->FillESD(digitsTree, clustersTree, esd);
1183         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1184       }
1185       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1186         fLoader[iDet]->UnloadRecPoints();
1187       }
1188
1189       if (fRawReader) {
1190         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1191       } else {
1192         reconstructor->FillESD(fRunLoader, esd);
1193       }
1194       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1195     }
1196   }
1197
1198   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1199     AliError(Form("the following detectors were not found: %s",
1200                   detStr.Data()));
1201     if (fStopOnError) return kFALSE;
1202   }
1203
1204   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1205                stopwatch.RealTime(),stopwatch.CpuTime()));
1206
1207   return kTRUE;
1208 }
1209
1210 //_____________________________________________________________________________
1211 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1212 {
1213   // Reads the trigger decision which is
1214   // stored in Trigger.root file and fills
1215   // the corresponding esd entries
1216
1217   AliInfo("Filling trigger information into the ESD");
1218
1219   if (fRawReader) {
1220     AliCTPRawStream input(fRawReader);
1221     if (!input.Next()) {
1222       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1223       return kFALSE;
1224     }
1225     esd->SetTriggerMask(input.GetClassMask());
1226     esd->SetTriggerCluster(input.GetClusterMask());
1227   }
1228   else {
1229     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1230     if (runloader) {
1231       if (!runloader->LoadTrigger()) {
1232         AliCentralTrigger *aCTP = runloader->GetTrigger();
1233         esd->SetTriggerMask(aCTP->GetClassMask());
1234         esd->SetTriggerCluster(aCTP->GetClusterMask());
1235       }
1236       else {
1237         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1238         return kFALSE;
1239       }
1240     }
1241     else {
1242       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1243       return kFALSE;
1244     }
1245   }
1246
1247   return kTRUE;
1248 }
1249
1250
1251
1252
1253
1254 //_____________________________________________________________________________
1255 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1256 {
1257   //
1258   // Filling information from RawReader Header
1259   //
1260
1261   AliInfo("Filling information from RawReader Header");
1262   esd->SetTimeStamp(0);
1263   esd->SetEventType(0);
1264   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1265   if (eventHeader){
1266     esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1267     esd->SetEventType((eventHeader->Get("Type")));
1268   }
1269
1270   return kTRUE;
1271 }
1272
1273
1274 //_____________________________________________________________________________
1275 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1276 {
1277 // check whether detName is contained in detectors
1278 // if yes, it is removed from detectors
1279
1280   // check if all detectors are selected
1281   if ((detectors.CompareTo("ALL") == 0) ||
1282       detectors.BeginsWith("ALL ") ||
1283       detectors.EndsWith(" ALL") ||
1284       detectors.Contains(" ALL ")) {
1285     detectors = "ALL";
1286     return kTRUE;
1287   }
1288
1289   // search for the given detector
1290   Bool_t result = kFALSE;
1291   if ((detectors.CompareTo(detName) == 0) ||
1292       detectors.BeginsWith(detName+" ") ||
1293       detectors.EndsWith(" "+detName) ||
1294       detectors.Contains(" "+detName+" ")) {
1295     detectors.ReplaceAll(detName, "");
1296     result = kTRUE;
1297   }
1298
1299   // clean up the detectors string
1300   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1301   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1302   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1303
1304   return result;
1305 }
1306
1307 //_____________________________________________________________________________
1308 Bool_t AliReconstruction::InitRunLoader()
1309 {
1310 // get or create the run loader
1311
1312   if (gAlice) delete gAlice;
1313   gAlice = NULL;
1314
1315   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1316     // load all base libraries to get the loader classes
1317     TString libs = gSystem->GetLibraries();
1318     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1319       TString detName = fgkDetectorName[iDet];
1320       if (detName == "HLT") continue;
1321       if (libs.Contains("lib" + detName + "base.so")) continue;
1322       gSystem->Load("lib" + detName + "base.so");
1323     }
1324     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1325     if (!fRunLoader) {
1326       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1327       CleanUp();
1328       return kFALSE;
1329     }
1330     fRunLoader->CdGAFile();
1331     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1332       if (fRunLoader->LoadgAlice() == 0) {
1333         gAlice = fRunLoader->GetAliRun();
1334         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1335       }
1336     }
1337     if (!gAlice && !fRawReader) {
1338       AliError(Form("no gAlice object found in file %s",
1339                     fGAliceFileName.Data()));
1340       CleanUp();
1341       return kFALSE;
1342     }
1343
1344   } else {               // galice.root does not exist
1345     if (!fRawReader) {
1346       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1347       CleanUp();
1348       return kFALSE;
1349     }
1350     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1351                                     AliConfig::GetDefaultEventFolderName(),
1352                                     "recreate");
1353     if (!fRunLoader) {
1354       AliError(Form("could not create run loader in file %s",
1355                     fGAliceFileName.Data()));
1356       CleanUp();
1357       return kFALSE;
1358     }
1359     fRunLoader->MakeTree("E");
1360     Int_t iEvent = 0;
1361     while (fRawReader->NextEvent()) {
1362       fRunLoader->SetEventNumber(iEvent);
1363       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1364                                      iEvent, iEvent);
1365       fRunLoader->MakeTree("H");
1366       fRunLoader->TreeE()->Fill();
1367       iEvent++;
1368     }
1369     fRawReader->RewindEvents();
1370     fRunLoader->WriteHeader("OVERWRITE");
1371     fRunLoader->CdGAFile();
1372     fRunLoader->Write(0, TObject::kOverwrite);
1373 //    AliTracker::SetFieldMap(???);
1374   }
1375
1376   return kTRUE;
1377 }
1378
1379 //_____________________________________________________________________________
1380 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1381 {
1382 // get the reconstructor object and the loader for a detector
1383
1384   if (fReconstructor[iDet]) return fReconstructor[iDet];
1385
1386   // load the reconstructor object
1387   TPluginManager* pluginManager = gROOT->GetPluginManager();
1388   TString detName = fgkDetectorName[iDet];
1389   TString recName = "Ali" + detName + "Reconstructor";
1390   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1391
1392   if (detName == "HLT") {
1393     if (!gROOT->GetClass("AliLevel3")) {
1394       gSystem->Load("libAliL3Src.so");
1395       gSystem->Load("libAliL3Misc.so");
1396       gSystem->Load("libAliL3Hough.so");
1397       gSystem->Load("libAliL3Comp.so");
1398     }
1399   }
1400
1401   AliReconstructor* reconstructor = NULL;
1402   // first check if a plugin is defined for the reconstructor
1403   TPluginHandler* pluginHandler =
1404     pluginManager->FindHandler("AliReconstructor", detName);
1405   // if not, add a plugin for it
1406   if (!pluginHandler) {
1407     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1408     TString libs = gSystem->GetLibraries();
1409     if (libs.Contains("lib" + detName + "base.so") ||
1410         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1411       pluginManager->AddHandler("AliReconstructor", detName,
1412                                 recName, detName + "rec", recName + "()");
1413     } else {
1414       pluginManager->AddHandler("AliReconstructor", detName,
1415                                 recName, detName, recName + "()");
1416     }
1417     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1418   }
1419   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1420     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1421   }
1422   if (reconstructor) {
1423     TObject* obj = fOptions.FindObject(detName.Data());
1424     if (obj) reconstructor->SetOption(obj->GetTitle());
1425     reconstructor->Init(fRunLoader);
1426     fReconstructor[iDet] = reconstructor;
1427   }
1428
1429   // get or create the loader
1430   if (detName != "HLT") {
1431     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1432     if (!fLoader[iDet]) {
1433       AliConfig::Instance()
1434         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1435                                 detName, detName);
1436       // first check if a plugin is defined for the loader
1437       TPluginHandler* pluginHandler =
1438         pluginManager->FindHandler("AliLoader", detName);
1439       // if not, add a plugin for it
1440       if (!pluginHandler) {
1441         TString loaderName = "Ali" + detName + "Loader";
1442         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1443         pluginManager->AddHandler("AliLoader", detName,
1444                                   loaderName, detName + "base",
1445                                   loaderName + "(const char*, TFolder*)");
1446         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1447       }
1448       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1449         fLoader[iDet] =
1450           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1451                                                  fRunLoader->GetEventFolder());
1452       }
1453       if (!fLoader[iDet]) {   // use default loader
1454         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1455       }
1456       if (!fLoader[iDet]) {
1457         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1458         if (fStopOnError) return NULL;
1459       } else {
1460         fRunLoader->AddLoader(fLoader[iDet]);
1461         fRunLoader->CdGAFile();
1462         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1463         fRunLoader->Write(0, TObject::kOverwrite);
1464       }
1465     }
1466   }
1467
1468   return reconstructor;
1469 }
1470
1471 //_____________________________________________________________________________
1472 Bool_t AliReconstruction::CreateVertexer()
1473 {
1474 // create the vertexer
1475
1476   fVertexer = NULL;
1477   AliReconstructor* itsReconstructor = GetReconstructor(0);
1478   if (itsReconstructor) {
1479     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1480   }
1481   if (!fVertexer) {
1482     AliWarning("couldn't create a vertexer for ITS");
1483     if (fStopOnError) return kFALSE;
1484   }
1485
1486   return kTRUE;
1487 }
1488
1489 //_____________________________________________________________________________
1490 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1491 {
1492 // create the trackers
1493
1494   TString detStr = detectors;
1495   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1496     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1497     AliReconstructor* reconstructor = GetReconstructor(iDet);
1498     if (!reconstructor) continue;
1499     TString detName = fgkDetectorName[iDet];
1500     if (detName == "HLT") {
1501       fRunHLTTracking = kTRUE;
1502       continue;
1503     }
1504
1505     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1506     if (!fTracker[iDet] && (iDet < 7)) {
1507       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1508       if (fStopOnError) return kFALSE;
1509     }
1510   }
1511
1512   return kTRUE;
1513 }
1514
1515 //_____________________________________________________________________________
1516 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1517 {
1518 // delete trackers and the run loader and close and delete the file
1519
1520   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1521     delete fReconstructor[iDet];
1522     fReconstructor[iDet] = NULL;
1523     fLoader[iDet] = NULL;
1524     delete fTracker[iDet];
1525     fTracker[iDet] = NULL;
1526   }
1527   delete fVertexer;
1528   fVertexer = NULL;
1529   delete fDiamondProfile;
1530   fDiamondProfile = NULL;
1531
1532   delete fRunLoader;
1533   fRunLoader = NULL;
1534   delete fRawReader;
1535   fRawReader = NULL;
1536
1537   if (file) {
1538     file->Close();
1539     delete file;
1540   }
1541
1542   if (fileOld) {
1543     fileOld->Close();
1544     delete fileOld;
1545     gSystem->Unlink("AliESDs.old.root");
1546   }
1547 }
1548
1549
1550 //_____________________________________________________________________________
1551 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1552 {
1553 // read the ESD event from a file
1554
1555   if (!esd) return kFALSE;
1556   char fileName[256];
1557   sprintf(fileName, "ESD_%d.%d_%s.root",
1558           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1559   if (gSystem->AccessPathName(fileName)) return kFALSE;
1560
1561   AliInfo(Form("reading ESD from file %s", fileName));
1562   AliDebug(1, Form("reading ESD from file %s", fileName));
1563   TFile* file = TFile::Open(fileName);
1564   if (!file || !file->IsOpen()) {
1565     AliError(Form("opening %s failed", fileName));
1566     delete file;
1567     return kFALSE;
1568   }
1569
1570   gROOT->cd();
1571   delete esd;
1572   esd = (AliESD*) file->Get("ESD");
1573   file->Close();
1574   delete file;
1575   return kTRUE;
1576 }
1577
1578 //_____________________________________________________________________________
1579 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1580 {
1581 // write the ESD event to a file
1582
1583   if (!esd) return;
1584   char fileName[256];
1585   sprintf(fileName, "ESD_%d.%d_%s.root",
1586           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1587
1588   AliDebug(1, Form("writing ESD to file %s", fileName));
1589   TFile* file = TFile::Open(fileName, "recreate");
1590   if (!file || !file->IsOpen()) {
1591     AliError(Form("opening %s failed", fileName));
1592   } else {
1593     esd->Write("ESD");
1594     file->Close();
1595   }
1596   delete file;
1597 }
1598
1599
1600
1601
1602 //_____________________________________________________________________________
1603 void AliReconstruction::CreateTag(TFile* file)
1604 {
1605   /////////////
1606   //muon code//
1607   ////////////
1608   Double_t fMUONMASS = 0.105658369;
1609   //Variables
1610   Double_t fX,fY,fZ ;
1611   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1612   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1613   Int_t fCharge;
1614   TLorentzVector fEPvector;
1615
1616   Float_t fZVertexCut = 10.0;
1617   Float_t fRhoVertexCut = 2.0;
1618
1619   Float_t fLowPtCut = 1.0;
1620   Float_t fHighPtCut = 3.0;
1621   Float_t fVeryHighPtCut = 10.0;
1622   ////////////
1623
1624   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1625
1626   // Creates the tags for all the events in a given ESD file
1627   Int_t ntrack;
1628   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1629   Int_t nPos, nNeg, nNeutr;
1630   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1631   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1632   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1633   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1634   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1635   Int_t fVertexflag;
1636   Int_t iRunNumber = 0;
1637   TString fVertexName("default");
1638
1639   AliRunTag *tag = new AliRunTag();
1640   AliEventTag *evTag = new AliEventTag();
1641   TTree ttag("T","A Tree with event tags");
1642   TBranch * btag = ttag.Branch("AliTAG", &tag);
1643   btag->SetCompressionLevel(9);
1644
1645   AliInfo(Form("Creating the tags......."));
1646
1647   if (!file || !file->IsOpen()) {
1648     AliError(Form("opening failed"));
1649     delete file;
1650     return ;
1651   }
1652   Int_t firstEvent = 0,lastEvent = 0;
1653   TTree *t = (TTree*) file->Get("esdTree");
1654   TBranch * b = t->GetBranch("ESD");
1655   AliESD *esd = 0;
1656   b->SetAddress(&esd);
1657
1658   b->GetEntry(0);
1659   Int_t iInitRunNumber = esd->GetRunNumber();
1660
1661   Int_t iNumberOfEvents = b->GetEntries();
1662   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1663     ntrack = 0;
1664     nPos = 0;
1665     nNeg = 0;
1666     nNeutr =0;
1667     nK0s = 0;
1668     nNeutrons = 0;
1669     nPi0s = 0;
1670     nGamas = 0;
1671     nProtons = 0;
1672     nKaons = 0;
1673     nPions = 0;
1674     nMuons = 0;
1675     nElectrons = 0;
1676     nCh1GeV = 0;
1677     nCh3GeV = 0;
1678     nCh10GeV = 0;
1679     nMu1GeV = 0;
1680     nMu3GeV = 0;
1681     nMu10GeV = 0;
1682     nEl1GeV = 0;
1683     nEl3GeV = 0;
1684     nEl10GeV = 0;
1685     maxPt = .0;
1686     meanPt = .0;
1687     totalP = .0;
1688     fVertexflag = 0;
1689
1690     b->GetEntry(iEventNumber);
1691     iRunNumber = esd->GetRunNumber();
1692     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1693     const AliESDVertex * vertexIn = esd->GetVertex();
1694     if (!vertexIn) AliError("ESD has not defined vertex.");
1695     if (vertexIn) fVertexName = vertexIn->GetName();
1696     if(fVertexName != "default") fVertexflag = 1;
1697     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1698       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1699       UInt_t status = esdTrack->GetStatus();
1700
1701       //select only tracks with ITS refit
1702       if ((status&AliESDtrack::kITSrefit)==0) continue;
1703       //select only tracks with TPC refit
1704       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1705
1706       //select only tracks with the "combined PID"
1707       if ((status&AliESDtrack::kESDpid)==0) continue;
1708       Double_t p[3];
1709       esdTrack->GetPxPyPz(p);
1710       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1711       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1712       totalP += momentum;
1713       meanPt += fPt;
1714       if(fPt > maxPt) maxPt = fPt;
1715
1716       if(esdTrack->GetSign() > 0) {
1717         nPos++;
1718         if(fPt > fLowPtCut) nCh1GeV++;
1719         if(fPt > fHighPtCut) nCh3GeV++;
1720         if(fPt > fVeryHighPtCut) nCh10GeV++;
1721       }
1722       if(esdTrack->GetSign() < 0) {
1723         nNeg++;
1724         if(fPt > fLowPtCut) nCh1GeV++;
1725         if(fPt > fHighPtCut) nCh3GeV++;
1726         if(fPt > fVeryHighPtCut) nCh10GeV++;
1727       }
1728       if(esdTrack->GetSign() == 0) nNeutr++;
1729
1730       //PID
1731       Double_t prob[5];
1732       esdTrack->GetESDpid(prob);
1733
1734       Double_t rcc = 0.0;
1735       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1736       if(rcc == 0.0) continue;
1737       //Bayes' formula
1738       Double_t w[5];
1739       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1740
1741       //protons
1742       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1743       //kaons
1744       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1745       //pions
1746       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1747       //electrons
1748       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1749         nElectrons++;
1750         if(fPt > fLowPtCut) nEl1GeV++;
1751         if(fPt > fHighPtCut) nEl3GeV++;
1752         if(fPt > fVeryHighPtCut) nEl10GeV++;
1753       }
1754       ntrack++;
1755     }//track loop
1756
1757     /////////////
1758     //muon code//
1759     ////////////
1760     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1761     // loop over all reconstructed tracks (also first track of combination)
1762     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1763       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1764       if (muonTrack == 0x0) continue;
1765
1766       // Coordinates at vertex
1767       fZ = muonTrack->GetZ();
1768       fY = muonTrack->GetBendingCoor();
1769       fX = muonTrack->GetNonBendingCoor();
1770
1771       fThetaX = muonTrack->GetThetaX();
1772       fThetaY = muonTrack->GetThetaY();
1773
1774       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1775       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1776       fPxRec = fPzRec * TMath::Tan(fThetaX);
1777       fPyRec = fPzRec * TMath::Tan(fThetaY);
1778       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1779
1780       //ChiSquare of the track if needed
1781       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1782       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1783       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1784
1785       // total number of muons inside a vertex cut
1786       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1787         nMuons++;
1788         if(fEPvector.Pt() > fLowPtCut) {
1789           nMu1GeV++;
1790           if(fEPvector.Pt() > fHighPtCut) {
1791             nMu3GeV++;
1792             if (fEPvector.Pt() > fVeryHighPtCut) {
1793               nMu10GeV++;
1794             }
1795           }
1796         }
1797       }
1798     }//muon track loop
1799
1800     // Fill the event tags
1801     if(ntrack != 0)
1802       meanPt = meanPt/ntrack;
1803
1804     evTag->SetEventId(iEventNumber+1);
1805     if (vertexIn) {
1806       evTag->SetVertexX(vertexIn->GetXv());
1807       evTag->SetVertexY(vertexIn->GetYv());
1808       evTag->SetVertexZ(vertexIn->GetZv());
1809       evTag->SetVertexZError(vertexIn->GetZRes());
1810     }
1811     evTag->SetVertexFlag(fVertexflag);
1812
1813     evTag->SetT0VertexZ(esd->GetT0zVertex());
1814
1815     evTag->SetTriggerMask(esd->GetTriggerMask());
1816     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1817
1818     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1819     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1820     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1821     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1822     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1823     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1824
1825
1826     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1827     evTag->SetNumOfPosTracks(nPos);
1828     evTag->SetNumOfNegTracks(nNeg);
1829     evTag->SetNumOfNeutrTracks(nNeutr);
1830
1831     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1832     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1833     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1834     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1835
1836     evTag->SetNumOfProtons(nProtons);
1837     evTag->SetNumOfKaons(nKaons);
1838     evTag->SetNumOfPions(nPions);
1839     evTag->SetNumOfMuons(nMuons);
1840     evTag->SetNumOfElectrons(nElectrons);
1841     evTag->SetNumOfPhotons(nGamas);
1842     evTag->SetNumOfPi0s(nPi0s);
1843     evTag->SetNumOfNeutrons(nNeutrons);
1844     evTag->SetNumOfKaon0s(nK0s);
1845
1846     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1847     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1848     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1849     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1850     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1851     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1852     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1853     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1854     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1855
1856     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1857     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1858
1859     evTag->SetTotalMomentum(totalP);
1860     evTag->SetMeanPt(meanPt);
1861     evTag->SetMaxPt(maxPt);
1862
1863     tag->SetRunId(iInitRunNumber);
1864     tag->AddEventTag(*evTag);
1865   }
1866   lastEvent = iNumberOfEvents;
1867
1868   ttag.Fill();
1869   tag->Clear();
1870
1871   char fileName[256];
1872   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1873           tag->GetRunId(),firstEvent,lastEvent );
1874   AliInfo(Form("writing tags to file %s", fileName));
1875   AliDebug(1, Form("writing tags to file %s", fileName));
1876
1877   TFile* ftag = TFile::Open(fileName, "recreate");
1878   ftag->cd();
1879   ttag.Write();
1880   ftag->Close();
1881   file->cd();
1882   delete tag;
1883   delete evTag;
1884 }
1885
1886 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1887 {
1888   // Write space-points which are then used in the alignment procedures
1889   // For the moment only ITS, TRD and TPC
1890
1891   // Load TOF clusters
1892   if (fTracker[3]){
1893     fLoader[3]->LoadRecPoints("read");
1894     TTree* tree = fLoader[3]->TreeR();
1895     if (!tree) {
1896       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1897       return;
1898     }
1899     fTracker[3]->LoadClusters(tree);
1900   }
1901   Int_t ntracks = esd->GetNumberOfTracks();
1902   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1903     {
1904       AliESDtrack *track = esd->GetTrack(itrack);
1905       Int_t nsp = 0;
1906       Int_t idx[200];
1907       for (Int_t iDet = 3; iDet >= 0; iDet--)
1908         nsp += track->GetNcls(iDet);
1909       if (nsp) {
1910         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1911         track->SetTrackPointArray(sp);
1912         Int_t isptrack = 0;
1913         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1914           AliTracker *tracker = fTracker[iDet];
1915           if (!tracker) continue;
1916           Int_t nspdet = track->GetNcls(iDet);
1917           if (nspdet <= 0) continue;
1918           track->GetClusters(iDet,idx);
1919           AliTrackPoint p;
1920           Int_t isp = 0;
1921           Int_t isp2 = 0;
1922           while (isp < nspdet) {
1923             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1924             const Int_t kNTPCmax = 159;
1925             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
1926             if (!isvalid) continue;
1927             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1928           }
1929         }
1930       }
1931     }
1932   if (fTracker[3]){
1933     fTracker[3]->UnloadClusters();
1934     fLoader[3]->UnloadRecPoints();
1935   }
1936 }