STEER/AliReconstruction.cxx

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