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