STEER/AliReconstruction.cxx

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