STEER/AliReconstruction.cxx

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