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