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