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