STEER/AliReconstruction.cxx

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