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