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