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