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