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