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   tree->SetBranchStatus("ESDfriend*",0);
 770   tree->Write();
 771   hlttree->Write();
 772
 773   // Create tags for the events in the ESD tree (the ESD tree is always present)
 774   // In case of empty events the tags will contain dummy values
 775   CreateTag(file);
 776   CleanUp(file, fileOld);
 777
 778   return kTRUE;
 779 }
 780
 781
 782 //_____________________________________________________________________________
 783 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 784 {
 785 // run the local reconstruction
 786
 787   TStopwatch stopwatch;
 788   stopwatch.Start();
 789
 790   TString detStr = detectors;
 791   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 792     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 793     AliReconstructor* reconstructor = GetReconstructor(iDet);
 794     if (!reconstructor) continue;
 795     if (reconstructor->HasLocalReconstruction()) continue;
 796
 797     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 798     TStopwatch stopwatchDet;
 799     stopwatchDet.Start();
 800     if (fRawReader) {
 801       fRawReader->RewindEvents();
 802       reconstructor->Reconstruct(fRunLoader, fRawReader);
 803     } else {
 804       reconstructor->Reconstruct(fRunLoader);
 805     }
 806     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 807                  fgkDetectorName[iDet],
 808                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 809   }
 810
 811   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 812     AliError(Form("the following detectors were not found: %s",
 813                   detStr.Data()));
 814     if (fStopOnError) return kFALSE;
 815   }
 816
 817   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 818                stopwatch.RealTime(),stopwatch.CpuTime()));
 819
 820   return kTRUE;
 821 }
 822
 823 //_____________________________________________________________________________
 824 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 825 {
 826 // run the local reconstruction
 827
 828   TStopwatch stopwatch;
 829   stopwatch.Start();
 830
 831   TString detStr = detectors;
 832   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 833     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 834     AliReconstructor* reconstructor = GetReconstructor(iDet);
 835     if (!reconstructor) continue;
 836     AliLoader* loader = fLoader[iDet];
 837
 838     // conversion of digits
 839     if (fRawReader && reconstructor->HasDigitConversion()) {
 840       AliInfo(Form("converting raw data digits into root objects for %s",
 841                    fgkDetectorName[iDet]));
 842       TStopwatch stopwatchDet;
 843       stopwatchDet.Start();
 844       loader->LoadDigits("update");
 845       loader->CleanDigits();
 846       loader->MakeDigitsContainer();
 847       TTree* digitsTree = loader->TreeD();
 848       reconstructor->ConvertDigits(fRawReader, digitsTree);
 849       loader->WriteDigits("OVERWRITE");
 850       loader->UnloadDigits();
 851       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 852                    fgkDetectorName[iDet],
 853                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 854     }
 855
 856     // local reconstruction
 857     if (!reconstructor->HasLocalReconstruction()) continue;
 858     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 859     TStopwatch stopwatchDet;
 860     stopwatchDet.Start();
 861     loader->LoadRecPoints("update");
 862     loader->CleanRecPoints();
 863     loader->MakeRecPointsContainer();
 864     TTree* clustersTree = loader->TreeR();
 865     if (fRawReader && !reconstructor->HasDigitConversion()) {
 866       reconstructor->Reconstruct(fRawReader, clustersTree);
 867     } else {
 868       loader->LoadDigits("read");
 869       TTree* digitsTree = loader->TreeD();
 870       if (!digitsTree) {
 871         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 872         if (fStopOnError) return kFALSE;
 873       } else {
 874         reconstructor->Reconstruct(digitsTree, clustersTree);
 875       }
 876       loader->UnloadDigits();
 877     }
 878     loader->WriteRecPoints("OVERWRITE");
 879     loader->UnloadRecPoints();
 880     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
 881                     fgkDetectorName[iDet],
 882                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 883   }
 884
 885   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 886     AliError(Form("the following detectors were not found: %s",
 887                   detStr.Data()));
 888     if (fStopOnError) return kFALSE;
 889   }
 890
 891   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 892                stopwatch.RealTime(),stopwatch.CpuTime()));
 893
 894   return kTRUE;
 895 }
 896
 897 //_____________________________________________________________________________
 898 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 899 {
 900 // run the barrel tracking
 901
 902   TStopwatch stopwatch;
 903   stopwatch.Start();
 904
 905   AliESDVertex* vertex = NULL;
 906   Double_t vtxPos[3] = {0, 0, 0};
 907   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 908   TArrayF mcVertex(3);
 909   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 910     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 911     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 912   }
 913
 914   if (fVertexer) {
 915     AliInfo("running the ITS vertex finder");
 916     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 917     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 918     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 919     if(!vertex){
 920       AliWarning("Vertex not found");
 921       vertex = new AliESDVertex();
 922       vertex->SetName("default");
 923     }
 924     else {
 925       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 926       vertex->SetName("reconstructed");
 927     }
 928
 929   } else {
 930     AliInfo("getting the primary vertex from MC");
 931     vertex = new AliESDVertex(vtxPos, vtxErr);
 932   }
 933
 934   if (vertex) {
 935     vertex->GetXYZ(vtxPos);
 936     vertex->GetSigmaXYZ(vtxErr);
 937   } else {
 938     AliWarning("no vertex reconstructed");
 939     vertex = new AliESDVertex(vtxPos, vtxErr);
 940   }
 941   esd->SetVertex(vertex);
 942   // if SPD multiplicity has been determined, it is stored in the ESD
 943   AliMultiplicity *mult= fVertexer->GetMultiplicity();
 944   if(mult)esd->SetMultiplicity(mult);
 945
 946   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 947     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
 948   }
 949   delete vertex;
 950
 951   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 952                stopwatch.RealTime(),stopwatch.CpuTime()));
 953
 954   return kTRUE;
 955 }
 956
 957 //_____________________________________________________________________________
 958 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
 959 {
 960 // run the HLT barrel tracking
 961
 962   TStopwatch stopwatch;
 963   stopwatch.Start();
 964
 965   if (!fRunLoader) {
 966     AliError("Missing runLoader!");
 967     return kFALSE;
 968   }
 969
 970   AliInfo("running HLT tracking");
 971
 972   // Get a pointer to the HLT reconstructor
 973   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
 974   if (!reconstructor) return kFALSE;
 975
 976   // TPC + ITS
 977   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 978     TString detName = fgkDetectorName[iDet];
 979     AliDebug(1, Form("%s HLT tracking", detName.Data()));
 980     reconstructor->SetOption(detName.Data());
 981     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
 982     if (!tracker) {
 983       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
 984       if (fStopOnError) return kFALSE;
 985       continue;
 986     }
 987     Double_t vtxPos[3];
 988     Double_t vtxErr[3]={0.005,0.005,0.010};
 989     const AliESDVertex *vertex = esd->GetVertex();
 990     vertex->GetXYZ(vtxPos);
 991     tracker->SetVertex(vtxPos,vtxErr);
 992     if(iDet != 1) {
 993       fLoader[iDet]->LoadRecPoints("read");
 994       TTree* tree = fLoader[iDet]->TreeR();
 995       if (!tree) {
 996         AliError(Form("Can't get the %s cluster tree", detName.Data()));
 997         return kFALSE;
 998       }
 999       tracker->LoadClusters(tree);
1000     }
1001     if (tracker->Clusters2Tracks(esd) != 0) {
1002       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1003       return kFALSE;
1004     }
1005     if(iDet != 1) {
1006       tracker->UnloadClusters();
1007     }
1008     delete tracker;
1009   }
1010
1011   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1012                stopwatch.RealTime(),stopwatch.CpuTime()));
1013
1014   return kTRUE;
1015 }
1016
1017 //_____________________________________________________________________________
1018 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1019 {
1020 // run the barrel tracking
1021
1022   TStopwatch stopwatch;
1023   stopwatch.Start();
1024
1025   AliInfo("running tracking");
1026
1027   // pass 1: TPC + ITS inwards
1028   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1029     if (!fTracker[iDet]) continue;
1030     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1031
1032     // load clusters
1033     fLoader[iDet]->LoadRecPoints("read");
1034     TTree* tree = fLoader[iDet]->TreeR();
1035     if (!tree) {
1036       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1037       return kFALSE;
1038     }
1039     fTracker[iDet]->LoadClusters(tree);
1040
1041     // run tracking
1042     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1043       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1044       return kFALSE;
1045     }
1046     if (fCheckPointLevel > 1) {
1047       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1048     }
1049     // preliminary PID in TPC needed by the ITS tracker
1050     if (iDet == 1) {
1051       GetReconstructor(1)->FillESD(fRunLoader, esd);
1052       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1053       AliESDpid::MakePID(esd);
1054     }
1055   }
1056
1057   // pass 2: ALL backwards
1058   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1059     if (!fTracker[iDet]) continue;
1060     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1061
1062     // load clusters
1063     if (iDet > 1) {     // all except ITS, TPC
1064       TTree* tree = NULL;
1065       fLoader[iDet]->LoadRecPoints("read");
1066       tree = fLoader[iDet]->TreeR();
1067       if (!tree) {
1068         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1069         return kFALSE;
1070       }
1071       fTracker[iDet]->LoadClusters(tree);
1072     }
1073
1074     // run tracking
1075     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1076       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1077       return kFALSE;
1078     }
1079     if (fCheckPointLevel > 1) {
1080       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1081     }
1082
1083     // unload clusters
1084     if (iDet > 2) {     // all except ITS, TPC, TRD
1085       fTracker[iDet]->UnloadClusters();
1086       fLoader[iDet]->UnloadRecPoints();
1087     }
1088     // updated PID in TPC needed by the ITS tracker -MI
1089     if (iDet == 1) {
1090       GetReconstructor(1)->FillESD(fRunLoader, esd);
1091       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1092       AliESDpid::MakePID(esd);
1093     }
1094   }
1095
1096   // write space-points to the ESD in case alignment data output
1097   // is switched on
1098   if (fWriteAlignmentData)
1099     WriteAlignmentData(esd);
1100
1101   // pass 3: TRD + TPC + ITS refit inwards
1102   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1103     if (!fTracker[iDet]) continue;
1104     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1105
1106     // run tracking
1107     if (fTracker[iDet]->RefitInward(esd) != 0) {
1108       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1109       return kFALSE;
1110     }
1111     if (fCheckPointLevel > 1) {
1112       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1113     }
1114
1115     // unload clusters
1116     fTracker[iDet]->UnloadClusters();
1117     fLoader[iDet]->UnloadRecPoints();
1118   }
1119   //
1120   // Propagate track to the vertex - if not done by ITS
1121   //
1122   Int_t ntracks = esd->GetNumberOfTracks();
1123   for (Int_t itrack=0; itrack<ntracks; itrack++){
1124     const Double_t kRadius  = 3;   // beam pipe radius
1125     const Double_t kMaxStep = 5;   // max step
1126     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1127     Double_t       fieldZ   = AliTracker::GetBz();  //
1128     AliESDtrack * track = esd->GetTrack(itrack);
1129     if (!track) continue;
1130     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1131     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1132     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1133   }
1134
1135   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1136                stopwatch.RealTime(),stopwatch.CpuTime()));
1137
1138   return kTRUE;
1139 }
1140
1141 //_____________________________________________________________________________
1142 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1143 {
1144 // fill the event summary data
1145
1146   TStopwatch stopwatch;
1147   stopwatch.Start();
1148   AliInfo("filling ESD");
1149
1150   TString detStr = detectors;
1151   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1152     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1153     AliReconstructor* reconstructor = GetReconstructor(iDet);
1154     if (!reconstructor) continue;
1155
1156     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1157       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1158       TTree* clustersTree = NULL;
1159       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1160         fLoader[iDet]->LoadRecPoints("read");
1161         clustersTree = fLoader[iDet]->TreeR();
1162         if (!clustersTree) {
1163           AliError(Form("Can't get the %s clusters tree",
1164                         fgkDetectorName[iDet]));
1165           if (fStopOnError) return kFALSE;
1166         }
1167       }
1168       if (fRawReader && !reconstructor->HasDigitConversion()) {
1169         reconstructor->FillESD(fRawReader, clustersTree, esd);
1170       } else {
1171         TTree* digitsTree = NULL;
1172         if (fLoader[iDet]) {
1173           fLoader[iDet]->LoadDigits("read");
1174           digitsTree = fLoader[iDet]->TreeD();
1175           if (!digitsTree) {
1176             AliError(Form("Can't get the %s digits tree",
1177                           fgkDetectorName[iDet]));
1178             if (fStopOnError) return kFALSE;
1179           }
1180         }
1181         reconstructor->FillESD(digitsTree, clustersTree, esd);
1182         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1183       }
1184       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1185         fLoader[iDet]->UnloadRecPoints();
1186       }
1187
1188       if (fRawReader) {
1189         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1190       } else {
1191         reconstructor->FillESD(fRunLoader, esd);
1192       }
1193       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1194     }
1195   }
1196
1197   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1198     AliError(Form("the following detectors were not found: %s",
1199                   detStr.Data()));
1200     if (fStopOnError) return kFALSE;
1201   }
1202
1203   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1204                stopwatch.RealTime(),stopwatch.CpuTime()));
1205
1206   return kTRUE;
1207 }
1208
1209 //_____________________________________________________________________________
1210 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1211 {
1212   // Reads the trigger decision which is
1213   // stored in Trigger.root file and fills
1214   // the corresponding esd entries
1215
1216   AliInfo("Filling trigger information into the ESD");
1217
1218   if (fRawReader) {
1219     AliCTPRawStream input(fRawReader);
1220     if (!input.Next()) {
1221       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1222       return kFALSE;
1223     }
1224     esd->SetTriggerMask(input.GetClassMask());
1225     esd->SetTriggerCluster(input.GetClusterMask());
1226   }
1227   else {
1228     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1229     if (runloader) {
1230       if (!runloader->LoadTrigger()) {
1231         AliCentralTrigger *aCTP = runloader->GetTrigger();
1232         esd->SetTriggerMask(aCTP->GetClassMask());
1233         esd->SetTriggerCluster(aCTP->GetClusterMask());
1234       }
1235       else {
1236         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1237         return kFALSE;
1238       }
1239     }
1240     else {
1241       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1242       return kFALSE;
1243     }
1244   }
1245
1246   return kTRUE;
1247 }
1248
1249
1250
1251
1252
1253 //_____________________________________________________________________________
1254 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1255 {
1256   //
1257   // Filling information from RawReader Header
1258   //
1259
1260   AliInfo("Filling information from RawReader Header");
1261   esd->SetTimeStamp(0);
1262   esd->SetEventType(0);
1263   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1264   if (eventHeader){
1265     esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1266     esd->SetEventType((eventHeader->Get("Type")));
1267   }
1268
1269   return kTRUE;
1270 }
1271
1272
1273 //_____________________________________________________________________________
1274 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1275 {
1276 // check whether detName is contained in detectors
1277 // if yes, it is removed from detectors
1278
1279   // check if all detectors are selected
1280   if ((detectors.CompareTo("ALL") == 0) ||
1281       detectors.BeginsWith("ALL ") ||
1282       detectors.EndsWith(" ALL") ||
1283       detectors.Contains(" ALL ")) {
1284     detectors = "ALL";
1285     return kTRUE;
1286   }
1287
1288   // search for the given detector
1289   Bool_t result = kFALSE;
1290   if ((detectors.CompareTo(detName) == 0) ||
1291       detectors.BeginsWith(detName+" ") ||
1292       detectors.EndsWith(" "+detName) ||
1293       detectors.Contains(" "+detName+" ")) {
1294     detectors.ReplaceAll(detName, "");
1295     result = kTRUE;
1296   }
1297
1298   // clean up the detectors string
1299   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1300   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1301   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1302
1303   return result;
1304 }
1305
1306 //_____________________________________________________________________________
1307 Bool_t AliReconstruction::InitRunLoader()
1308 {
1309 // get or create the run loader
1310
1311   if (gAlice) delete gAlice;
1312   gAlice = NULL;
1313
1314   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1315     // load all base libraries to get the loader classes
1316     TString libs = gSystem->GetLibraries();
1317     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1318       TString detName = fgkDetectorName[iDet];
1319       if (detName == "HLT") continue;
1320       if (libs.Contains("lib" + detName + "base.so")) continue;
1321       gSystem->Load("lib" + detName + "base.so");
1322     }
1323     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1324     if (!fRunLoader) {
1325       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1326       CleanUp();
1327       return kFALSE;
1328     }
1329     fRunLoader->CdGAFile();
1330     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1331       if (fRunLoader->LoadgAlice() == 0) {
1332         gAlice = fRunLoader->GetAliRun();
1333         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1334       }
1335     }
1336     if (!gAlice && !fRawReader) {
1337       AliError(Form("no gAlice object found in file %s",
1338                     fGAliceFileName.Data()));
1339       CleanUp();
1340       return kFALSE;
1341     }
1342
1343   } else {               // galice.root does not exist
1344     if (!fRawReader) {
1345       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1346       CleanUp();
1347       return kFALSE;
1348     }
1349     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1350                                     AliConfig::GetDefaultEventFolderName(),
1351                                     "recreate");
1352     if (!fRunLoader) {
1353       AliError(Form("could not create run loader in file %s",
1354                     fGAliceFileName.Data()));
1355       CleanUp();
1356       return kFALSE;
1357     }
1358     fRunLoader->MakeTree("E");
1359     Int_t iEvent = 0;
1360     while (fRawReader->NextEvent()) {
1361       fRunLoader->SetEventNumber(iEvent);
1362       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1363                                      iEvent, iEvent);
1364       fRunLoader->MakeTree("H");
1365       fRunLoader->TreeE()->Fill();
1366       iEvent++;
1367     }
1368     fRawReader->RewindEvents();
1369     fRunLoader->WriteHeader("OVERWRITE");
1370     fRunLoader->CdGAFile();
1371     fRunLoader->Write(0, TObject::kOverwrite);
1372 //    AliTracker::SetFieldMap(???);
1373   }
1374
1375   return kTRUE;
1376 }
1377
1378 //_____________________________________________________________________________
1379 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1380 {
1381 // get the reconstructor object and the loader for a detector
1382
1383   if (fReconstructor[iDet]) return fReconstructor[iDet];
1384
1385   // load the reconstructor object
1386   TPluginManager* pluginManager = gROOT->GetPluginManager();
1387   TString detName = fgkDetectorName[iDet];
1388   TString recName = "Ali" + detName + "Reconstructor";
1389   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1390
1391   if (detName == "HLT") {
1392     if (!gROOT->GetClass("AliLevel3")) {
1393       gSystem->Load("libAliL3Src.so");
1394       gSystem->Load("libAliL3Misc.so");
1395       gSystem->Load("libAliL3Hough.so");
1396       gSystem->Load("libAliL3Comp.so");
1397     }
1398   }
1399
1400   AliReconstructor* reconstructor = NULL;
1401   // first check if a plugin is defined for the reconstructor
1402   TPluginHandler* pluginHandler =
1403     pluginManager->FindHandler("AliReconstructor", detName);
1404   // if not, add a plugin for it
1405   if (!pluginHandler) {
1406     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1407     TString libs = gSystem->GetLibraries();
1408     if (libs.Contains("lib" + detName + "base.so") ||
1409         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1410       pluginManager->AddHandler("AliReconstructor", detName,
1411                                 recName, detName + "rec", recName + "()");
1412     } else {
1413       pluginManager->AddHandler("AliReconstructor", detName,
1414                                 recName, detName, recName + "()");
1415     }
1416     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1417   }
1418   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1419     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1420   }
1421   if (reconstructor) {
1422     TObject* obj = fOptions.FindObject(detName.Data());
1423     if (obj) reconstructor->SetOption(obj->GetTitle());
1424     reconstructor->Init(fRunLoader);
1425     fReconstructor[iDet] = reconstructor;
1426   }
1427
1428   // get or create the loader
1429   if (detName != "HLT") {
1430     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1431     if (!fLoader[iDet]) {
1432       AliConfig::Instance()
1433         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1434                                 detName, detName);
1435       // first check if a plugin is defined for the loader
1436       TPluginHandler* pluginHandler =
1437         pluginManager->FindHandler("AliLoader", detName);
1438       // if not, add a plugin for it
1439       if (!pluginHandler) {
1440         TString loaderName = "Ali" + detName + "Loader";
1441         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1442         pluginManager->AddHandler("AliLoader", detName,
1443                                   loaderName, detName + "base",
1444                                   loaderName + "(const char*, TFolder*)");
1445         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1446       }
1447       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1448         fLoader[iDet] =
1449           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1450                                                  fRunLoader->GetEventFolder());
1451       }
1452       if (!fLoader[iDet]) {   // use default loader
1453         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1454       }
1455       if (!fLoader[iDet]) {
1456         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1457         if (fStopOnError) return NULL;
1458       } else {
1459         fRunLoader->AddLoader(fLoader[iDet]);
1460         fRunLoader->CdGAFile();
1461         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1462         fRunLoader->Write(0, TObject::kOverwrite);
1463       }
1464     }
1465   }
1466
1467   return reconstructor;
1468 }
1469
1470 //_____________________________________________________________________________
1471 Bool_t AliReconstruction::CreateVertexer()
1472 {
1473 // create the vertexer
1474
1475   fVertexer = NULL;
1476   AliReconstructor* itsReconstructor = GetReconstructor(0);
1477   if (itsReconstructor) {
1478     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1479   }
1480   if (!fVertexer) {
1481     AliWarning("couldn't create a vertexer for ITS");
1482     if (fStopOnError) return kFALSE;
1483   }
1484
1485   return kTRUE;
1486 }
1487
1488 //_____________________________________________________________________________
1489 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1490 {
1491 // create the trackers
1492
1493   TString detStr = detectors;
1494   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1495     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1496     AliReconstructor* reconstructor = GetReconstructor(iDet);
1497     if (!reconstructor) continue;
1498     TString detName = fgkDetectorName[iDet];
1499     if (detName == "HLT") {
1500       fRunHLTTracking = kTRUE;
1501       continue;
1502     }
1503
1504     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1505     if (!fTracker[iDet] && (iDet < 7)) {
1506       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1507       if (fStopOnError) return kFALSE;
1508     }
1509   }
1510
1511   return kTRUE;
1512 }
1513
1514 //_____________________________________________________________________________
1515 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1516 {
1517 // delete trackers and the run loader and close and delete the file
1518
1519   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1520     delete fReconstructor[iDet];
1521     fReconstructor[iDet] = NULL;
1522     fLoader[iDet] = NULL;
1523     delete fTracker[iDet];
1524     fTracker[iDet] = NULL;
1525   }
1526   delete fVertexer;
1527   fVertexer = NULL;
1528   delete fDiamondProfile;
1529   fDiamondProfile = NULL;
1530
1531   delete fRunLoader;
1532   fRunLoader = NULL;
1533   delete fRawReader;
1534   fRawReader = NULL;
1535
1536   if (file) {
1537     file->Close();
1538     delete file;
1539   }
1540
1541   if (fileOld) {
1542     fileOld->Close();
1543     delete fileOld;
1544     gSystem->Unlink("AliESDs.old.root");
1545   }
1546 }
1547
1548
1549 //_____________________________________________________________________________
1550 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1551 {
1552 // read the ESD event from a file
1553
1554   if (!esd) return kFALSE;
1555   char fileName[256];
1556   sprintf(fileName, "ESD_%d.%d_%s.root",
1557           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1558   if (gSystem->AccessPathName(fileName)) return kFALSE;
1559
1560   AliInfo(Form("reading ESD from file %s", fileName));
1561   AliDebug(1, Form("reading ESD from file %s", fileName));
1562   TFile* file = TFile::Open(fileName);
1563   if (!file || !file->IsOpen()) {
1564     AliError(Form("opening %s failed", fileName));
1565     delete file;
1566     return kFALSE;
1567   }
1568
1569   gROOT->cd();
1570   delete esd;
1571   esd = (AliESD*) file->Get("ESD");
1572   file->Close();
1573   delete file;
1574   return kTRUE;
1575 }
1576
1577 //_____________________________________________________________________________
1578 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1579 {
1580 // write the ESD event to a file
1581
1582   if (!esd) return;
1583   char fileName[256];
1584   sprintf(fileName, "ESD_%d.%d_%s.root",
1585           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1586
1587   AliDebug(1, Form("writing ESD to file %s", fileName));
1588   TFile* file = TFile::Open(fileName, "recreate");
1589   if (!file || !file->IsOpen()) {
1590     AliError(Form("opening %s failed", fileName));
1591   } else {
1592     esd->Write("ESD");
1593     file->Close();
1594   }
1595   delete file;
1596 }
1597
1598
1599
1600
1601 //_____________________________________________________________________________
1602 void AliReconstruction::CreateTag(TFile* file)
1603 {
1604   /////////////
1605   //muon code//
1606   ////////////
1607   Double_t fMUONMASS = 0.105658369;
1608   //Variables
1609   Double_t fX,fY,fZ ;
1610   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1611   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1612   Int_t fCharge;
1613   TLorentzVector fEPvector;
1614
1615   Float_t fZVertexCut = 10.0;
1616   Float_t fRhoVertexCut = 2.0;
1617
1618   Float_t fLowPtCut = 1.0;
1619   Float_t fHighPtCut = 3.0;
1620   Float_t fVeryHighPtCut = 10.0;
1621   ////////////
1622
1623   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1624
1625   // Creates the tags for all the events in a given ESD file
1626   Int_t ntrack;
1627   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1628   Int_t nPos, nNeg, nNeutr;
1629   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1630   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1631   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1632   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1633   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1634   Int_t fVertexflag;
1635   Int_t iRunNumber = 0;
1636   TString fVertexName("default");
1637
1638   AliRunTag *tag = new AliRunTag();
1639   AliEventTag *evTag = new AliEventTag();
1640   TTree ttag("T","A Tree with event tags");
1641   TBranch * btag = ttag.Branch("AliTAG", &tag);
1642   btag->SetCompressionLevel(9);
1643
1644   AliInfo(Form("Creating the tags......."));
1645
1646   if (!file || !file->IsOpen()) {
1647     AliError(Form("opening failed"));
1648     delete file;
1649     return ;
1650   }
1651   Int_t firstEvent = 0,lastEvent = 0;
1652   TTree *t = (TTree*) file->Get("esdTree");
1653   TBranch * b = t->GetBranch("ESD");
1654   AliESD *esd = 0;
1655   b->SetAddress(&esd);
1656
1657   b->GetEntry(0);
1658   Int_t iInitRunNumber = esd->GetRunNumber();
1659
1660   Int_t iNumberOfEvents = b->GetEntries();
1661   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1662     ntrack = 0;
1663     nPos = 0;
1664     nNeg = 0;
1665     nNeutr =0;
1666     nK0s = 0;
1667     nNeutrons = 0;
1668     nPi0s = 0;
1669     nGamas = 0;
1670     nProtons = 0;
1671     nKaons = 0;
1672     nPions = 0;
1673     nMuons = 0;
1674     nElectrons = 0;
1675     nCh1GeV = 0;
1676     nCh3GeV = 0;
1677     nCh10GeV = 0;
1678     nMu1GeV = 0;
1679     nMu3GeV = 0;
1680     nMu10GeV = 0;
1681     nEl1GeV = 0;
1682     nEl3GeV = 0;
1683     nEl10GeV = 0;
1684     maxPt = .0;
1685     meanPt = .0;
1686     totalP = .0;
1687     fVertexflag = 0;
1688
1689     b->GetEntry(iEventNumber);
1690     iRunNumber = esd->GetRunNumber();
1691     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1692     const AliESDVertex * vertexIn = esd->GetVertex();
1693     if (!vertexIn) AliError("ESD has not defined vertex.");
1694     if (vertexIn) fVertexName = vertexIn->GetName();
1695     if(fVertexName != "default") fVertexflag = 1;
1696     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1697       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1698       UInt_t status = esdTrack->GetStatus();
1699
1700       //select only tracks with ITS refit
1701       if ((status&AliESDtrack::kITSrefit)==0) continue;
1702       //select only tracks with TPC refit
1703       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1704
1705       //select only tracks with the "combined PID"
1706       if ((status&AliESDtrack::kESDpid)==0) continue;
1707       Double_t p[3];
1708       esdTrack->GetPxPyPz(p);
1709       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1710       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1711       totalP += momentum;
1712       meanPt += fPt;
1713       if(fPt > maxPt) maxPt = fPt;
1714
1715       if(esdTrack->GetSign() > 0) {
1716         nPos++;
1717         if(fPt > fLowPtCut) nCh1GeV++;
1718         if(fPt > fHighPtCut) nCh3GeV++;
1719         if(fPt > fVeryHighPtCut) nCh10GeV++;
1720       }
1721       if(esdTrack->GetSign() < 0) {
1722         nNeg++;
1723         if(fPt > fLowPtCut) nCh1GeV++;
1724         if(fPt > fHighPtCut) nCh3GeV++;
1725         if(fPt > fVeryHighPtCut) nCh10GeV++;
1726       }
1727       if(esdTrack->GetSign() == 0) nNeutr++;
1728
1729       //PID
1730       Double_t prob[5];
1731       esdTrack->GetESDpid(prob);
1732
1733       Double_t rcc = 0.0;
1734       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1735       if(rcc == 0.0) continue;
1736       //Bayes' formula
1737       Double_t w[5];
1738       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1739
1740       //protons
1741       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1742       //kaons
1743       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1744       //pions
1745       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1746       //electrons
1747       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1748         nElectrons++;
1749         if(fPt > fLowPtCut) nEl1GeV++;
1750         if(fPt > fHighPtCut) nEl3GeV++;
1751         if(fPt > fVeryHighPtCut) nEl10GeV++;
1752       }
1753       ntrack++;
1754     }//track loop
1755
1756     /////////////
1757     //muon code//
1758     ////////////
1759     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1760     // loop over all reconstructed tracks (also first track of combination)
1761     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1762       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1763       if (muonTrack == 0x0) continue;
1764
1765       // Coordinates at vertex
1766       fZ = muonTrack->GetZ();
1767       fY = muonTrack->GetBendingCoor();
1768       fX = muonTrack->GetNonBendingCoor();
1769
1770       fThetaX = muonTrack->GetThetaX();
1771       fThetaY = muonTrack->GetThetaY();
1772
1773       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1774       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1775       fPxRec = fPzRec * TMath::Tan(fThetaX);
1776       fPyRec = fPzRec * TMath::Tan(fThetaY);
1777       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1778
1779       //ChiSquare of the track if needed
1780       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1781       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1782       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1783
1784       // total number of muons inside a vertex cut
1785       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1786         nMuons++;
1787         if(fEPvector.Pt() > fLowPtCut) {
1788           nMu1GeV++;
1789           if(fEPvector.Pt() > fHighPtCut) {
1790             nMu3GeV++;
1791             if (fEPvector.Pt() > fVeryHighPtCut) {
1792               nMu10GeV++;
1793             }
1794           }
1795         }
1796       }
1797     }//muon track loop
1798
1799     // Fill the event tags
1800     if(ntrack != 0)
1801       meanPt = meanPt/ntrack;
1802
1803     evTag->SetEventId(iEventNumber+1);
1804     if (vertexIn) {
1805       evTag->SetVertexX(vertexIn->GetXv());
1806       evTag->SetVertexY(vertexIn->GetYv());
1807       evTag->SetVertexZ(vertexIn->GetZv());
1808       evTag->SetVertexZError(vertexIn->GetZRes());
1809     }
1810     evTag->SetVertexFlag(fVertexflag);
1811
1812     evTag->SetT0VertexZ(esd->GetT0zVertex());
1813
1814     evTag->SetTriggerMask(esd->GetTriggerMask());
1815     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1816
1817     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1818     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1819     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1820     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1821     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1822     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1823
1824
1825     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1826     evTag->SetNumOfPosTracks(nPos);
1827     evTag->SetNumOfNegTracks(nNeg);
1828     evTag->SetNumOfNeutrTracks(nNeutr);
1829
1830     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1831     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1832     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1833     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1834
1835     evTag->SetNumOfProtons(nProtons);
1836     evTag->SetNumOfKaons(nKaons);
1837     evTag->SetNumOfPions(nPions);
1838     evTag->SetNumOfMuons(nMuons);
1839     evTag->SetNumOfElectrons(nElectrons);
1840     evTag->SetNumOfPhotons(nGamas);
1841     evTag->SetNumOfPi0s(nPi0s);
1842     evTag->SetNumOfNeutrons(nNeutrons);
1843     evTag->SetNumOfKaon0s(nK0s);
1844
1845     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1846     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1847     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1848     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1849     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1850     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1851     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1852     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1853     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1854
1855     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1856     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1857
1858     evTag->SetTotalMomentum(totalP);
1859     evTag->SetMeanPt(meanPt);
1860     evTag->SetMaxPt(maxPt);
1861
1862     tag->SetRunId(iInitRunNumber);
1863     tag->AddEventTag(*evTag);
1864   }
1865   lastEvent = iNumberOfEvents;
1866
1867   ttag.Fill();
1868   tag->Clear();
1869
1870   char fileName[256];
1871   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1872           tag->GetRunId(),firstEvent,lastEvent );
1873   AliInfo(Form("writing tags to file %s", fileName));
1874   AliDebug(1, Form("writing tags to file %s", fileName));
1875
1876   TFile* ftag = TFile::Open(fileName, "recreate");
1877   ftag->cd();
1878   ttag.Write();
1879   ftag->Close();
1880   file->cd();
1881   delete tag;
1882   delete evTag;
1883 }
1884
1885 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1886 {
1887   // Write space-points which are then used in the alignment procedures
1888   // For the moment only ITS, TRD and TPC
1889
1890   // Load TOF clusters
1891   if (fTracker[3]){
1892     fLoader[3]->LoadRecPoints("read");
1893     TTree* tree = fLoader[3]->TreeR();
1894     if (!tree) {
1895       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1896       return;
1897     }
1898     fTracker[3]->LoadClusters(tree);
1899   }
1900   Int_t ntracks = esd->GetNumberOfTracks();
1901   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1902     {
1903       AliESDtrack *track = esd->GetTrack(itrack);
1904       Int_t nsp = 0;
1905       Int_t idx[200];
1906       for (Int_t iDet = 3; iDet >= 0; iDet--)
1907         nsp += track->GetNcls(iDet);
1908       if (nsp) {
1909         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1910         track->SetTrackPointArray(sp);
1911         Int_t isptrack = 0;
1912         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1913           AliTracker *tracker = fTracker[iDet];
1914           if (!tracker) continue;
1915           Int_t nspdet = track->GetNcls(iDet);
1916           if (nspdet <= 0) continue;
1917           track->GetClusters(iDet,idx);
1918           AliTrackPoint p;
1919           Int_t isp = 0;
1920           Int_t isp2 = 0;
1921           while (isp < nspdet) {
1922             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1923             const Int_t kNTPCmax = 159;
1924             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
1925             if (!isvalid) continue;
1926             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1927           }
1928         }
1929       }
1930     }
1931   if (fTracker[3]){
1932     fTracker[3]->UnloadClusters();
1933     fLoader[3]->UnloadRecPoints();
1934   }
1935 }