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