STEER/AliReconstruction.cxx

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