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