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