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
 143 ClassImp(AliReconstruction)
 144
 145
 146 //_____________________________________________________________________________
 147 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
 148
 149 //_____________________________________________________________________________
 150 AliReconstruction::AliReconstruction(const char* gAliceFilename,
 151                                      const char* name, const char* title) :
 152   TNamed(name, title),
 153
 154   fRunLocalReconstruction("ALL"),
 155   fUniformField(kTRUE),
 156   fRunVertexFinder(kTRUE),
 157   fRunHLTTracking(kFALSE),
 158   fRunTracking("ALL"),
 159   fFillESD("ALL"),
 160   fGAliceFileName(gAliceFilename),
 161   fInput(""),
 162   fFirstEvent(0),
 163   fLastEvent(-1),
 164   fStopOnError(kFALSE),
 165   fCheckPointLevel(0),
 166   fOptions(),
 167
 168   fRunLoader(NULL),
 169   fRawReader(NULL),
 170
 171   fVertexer(NULL),
 172
 173   fWriteAlignmentData(kFALSE)
 174 {
 175 // create reconstruction object with default parameters
 176
 177   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 178     fReconstructor[iDet] = NULL;
 179     fLoader[iDet] = NULL;
 180     fTracker[iDet] = NULL;
 181   }
 182   AliPID pid;
 183   // Import TGeo geometry
 184   TString geom(gSystem->DirName(gAliceFilename));
 185   geom += "/geometry.root";
 186   TGeoManager::Import(geom.Data());
 187 }
 188
 189 //_____________________________________________________________________________
 190 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
 191   TNamed(rec),
 192
 193   fRunLocalReconstruction(rec.fRunLocalReconstruction),
 194   fUniformField(rec.fUniformField),
 195   fRunVertexFinder(rec.fRunVertexFinder),
 196   fRunHLTTracking(rec.fRunHLTTracking),
 197   fRunTracking(rec.fRunTracking),
 198   fFillESD(rec.fFillESD),
 199   fGAliceFileName(rec.fGAliceFileName),
 200   fInput(rec.fInput),
 201   fFirstEvent(rec.fFirstEvent),
 202   fLastEvent(rec.fLastEvent),
 203   fStopOnError(rec.fStopOnError),
 204   fCheckPointLevel(0),
 205   fOptions(),
 206
 207   fRunLoader(NULL),
 208   fRawReader(NULL),
 209
 210   fVertexer(NULL),
 211
 212   fWriteAlignmentData(rec.fWriteAlignmentData)
 213 {
 214 // copy constructor
 215
 216   for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
 217     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
 218   }
 219   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 220     fReconstructor[iDet] = NULL;
 221     fLoader[iDet] = NULL;
 222     fTracker[iDet] = NULL;
 223   }
 224 }
 225
 226 //_____________________________________________________________________________
 227 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 228 {
 229 // assignment operator
 230
 231   this->~AliReconstruction();
 232   new(this) AliReconstruction(rec);
 233   return *this;
 234 }
 235
 236 //_____________________________________________________________________________
 237 AliReconstruction::~AliReconstruction()
 238 {
 239 // clean up
 240
 241   CleanUp();
 242   fOptions.Delete();
 243 }
 244
 245
 246 //_____________________________________________________________________________
 247 void AliReconstruction::SetGAliceFile(const char* fileName)
 248 {
 249 // set the name of the galice file
 250
 251   fGAliceFileName = fileName;
 252 }
 253
 254 //_____________________________________________________________________________
 255 void AliReconstruction::SetOption(const char* detector, const char* option)
 256 {
 257 // set options for the reconstruction of a detector
 258
 259   TObject* obj = fOptions.FindObject(detector);
 260   if (obj) fOptions.Remove(obj);
 261   fOptions.Add(new TNamed(detector, option));
 262 }
 263
 264
 265 //_____________________________________________________________________________
 266 Bool_t AliReconstruction::Run(const char* input,
 267                               Int_t firstEvent, Int_t lastEvent)
 268 {
 269 // run the reconstruction
 270
 271   // First check if we have any CDB storage set, because it is used
 272   // to retrieve the calibration and alignment constants
 273
 274   AliCDBManager* man = AliCDBManager::Instance();
 275   if (!man->IsDefaultStorageSet())
 276   {
 277     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 278     AliWarning("No default CDB storage set, so I will use $ALICE_ROOT");
 279     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 280     man->SetDefaultStorage("local://$ALICE_ROOT");
 281   }
 282
 283   // set the input
 284   if (!input) input = fInput.Data();
 285   TString fileName(input);
 286   if (fileName.EndsWith("/")) {
 287     fRawReader = new AliRawReaderFile(fileName);
 288   } else if (fileName.EndsWith(".root")) {
 289     fRawReader = new AliRawReaderRoot(fileName);
 290   } else if (!fileName.IsNull()) {
 291     fRawReader = new AliRawReaderDate(fileName);
 292     fRawReader->SelectEvents(7);
 293   }
 294
 295   // get the run loader
 296   if (!InitRunLoader()) return kFALSE;
 297
 298   // local reconstruction
 299   if (!fRunLocalReconstruction.IsNull()) {
 300     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
 301       if (fStopOnError) {CleanUp(); return kFALSE;}
 302     }
 303   }
 304 //  if (!fRunVertexFinder && fRunTracking.IsNull() &&
 305 //      fFillESD.IsNull()) return kTRUE;
 306
 307   // get vertexer
 308   if (fRunVertexFinder && !CreateVertexer()) {
 309     if (fStopOnError) {
 310       CleanUp();
 311       return kFALSE;
 312     }
 313   }
 314
 315   // get trackers
 316   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
 317     if (fStopOnError) {
 318       CleanUp();
 319       return kFALSE;
 320     }
 321   }
 322
 323
 324   TStopwatch stopwatch;
 325   stopwatch.Start();
 326
 327   // get the possibly already existing ESD file and tree
 328   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
 329   TFile* fileOld = NULL;
 330   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
 331   if (!gSystem->AccessPathName("AliESDs.root")){
 332     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
 333     fileOld = TFile::Open("AliESDs.old.root");
 334     if (fileOld && fileOld->IsOpen()) {
 335       treeOld = (TTree*) fileOld->Get("esdTree");
 336       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
 337       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
 338       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
 339     }
 340   }
 341
 342   // create the ESD output file and tree
 343   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 344   if (!file->IsOpen()) {
 345     AliError("opening AliESDs.root failed");
 346     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 347   }
 348   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 349   tree->Branch("ESD", "AliESD", &esd);
 350   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
 351   hlttree->Branch("ESD", "AliESD", &hltesd);
 352   delete esd; delete hltesd;
 353   esd = NULL; hltesd = NULL;
 354   gROOT->cd();
 355
 356   // loop over events
 357   if (fRawReader) fRawReader->RewindEvents();
 358
 359   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 360     if (fRawReader) fRawReader->NextEvent();
 361     if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
 362       // copy old ESD to the new one
 363       if (treeOld) {
 364         treeOld->SetBranchAddress("ESD", &esd);
 365         treeOld->GetEntry(iEvent);
 366       }
 367       tree->Fill();
 368       if (hlttreeOld) {
 369         hlttreeOld->SetBranchAddress("ESD", &hltesd);
 370         hlttreeOld->GetEntry(iEvent);
 371       }
 372       hlttree->Fill();
 373       continue;
 374     }
 375
 376     AliInfo(Form("processing event %d", iEvent));
 377     fRunLoader->GetEvent(iEvent);
 378
 379     char fileName[256];
 380     sprintf(fileName, "ESD_%d.%d_final.root",
 381             fRunLoader->GetHeader()->GetRun(),
 382             fRunLoader->GetHeader()->GetEventNrInRun());
 383     if (!gSystem->AccessPathName(fileName)) continue;
 384
 385     // local reconstruction
 386     if (!fRunLocalReconstruction.IsNull()) {
 387       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
 388         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 389       }
 390     }
 391
 392     esd = new AliESD; hltesd = new AliESD;
 393     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 394     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 395     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 396     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 397     if (gAlice) {
 398       esd->SetMagneticField(AliTracker::GetBz());
 399       hltesd->SetMagneticField(AliTracker::GetBz());
 400     } else {
 401       // ???
 402     }
 403
 404     // vertex finder
 405     if (fRunVertexFinder) {
 406       if (!ReadESD(esd, "vertex")) {
 407         if (!RunVertexFinder(esd)) {
 408           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 409         }
 410         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 411       }
 412     }
 413
 414     // HLT tracking
 415     if (!fRunTracking.IsNull()) {
 416       if (fRunHLTTracking) {
 417         hltesd->SetVertex(esd->GetVertex());
 418         if (!RunHLTTracking(hltesd)) {
 419           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 420         }
 421       }
 422     }
 423
 424     // barrel tracking
 425     if (!fRunTracking.IsNull()) {
 426       if (!ReadESD(esd, "tracking")) {
 427         if (!RunTracking(esd)) {
 428           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 429         }
 430         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 431       }
 432     }
 433
 434     // fill ESD
 435     if (!fFillESD.IsNull()) {
 436       if (!FillESD(esd, fFillESD)) {
 437         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 438       }
 439     }
 440
 441     // combined PID
 442     AliESDpid::MakePID(esd);
 443     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 444
 445     // write ESD
 446     tree->Fill();
 447     // write HLT ESD
 448     hlttree->Fill();
 449
 450     if (fCheckPointLevel > 0)  WriteESD(esd, "final");
 451
 452     delete esd; delete hltesd;
 453     esd = NULL; hltesd = NULL;
 454   }
 455
 456   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
 457                stopwatch.RealTime(),stopwatch.CpuTime()));
 458
 459   file->cd();
 460   tree->Write();
 461   hlttree->Write();
 462
 463   // Create tags for the events in the ESD tree (the ESD tree is always present)
 464   // In case of empty events the tags will contain dummy values
 465   CreateTag(file);
 466   CleanUp(file, fileOld);
 467
 468   return kTRUE;
 469 }
 470
 471
 472 //_____________________________________________________________________________
 473 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 474 {
 475 // run the local reconstruction
 476
 477   TStopwatch stopwatch;
 478   stopwatch.Start();
 479
 480   TString detStr = detectors;
 481   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 482     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 483     AliReconstructor* reconstructor = GetReconstructor(iDet);
 484     if (!reconstructor) continue;
 485     if (reconstructor->HasLocalReconstruction()) continue;
 486
 487     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 488     TStopwatch stopwatchDet;
 489     stopwatchDet.Start();
 490     if (fRawReader) {
 491       fRawReader->RewindEvents();
 492       reconstructor->Reconstruct(fRunLoader, fRawReader);
 493     } else {
 494       reconstructor->Reconstruct(fRunLoader);
 495     }
 496     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 497                  fgkDetectorName[iDet],
 498                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 499   }
 500
 501   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 502     AliError(Form("the following detectors were not found: %s",
 503                   detStr.Data()));
 504     if (fStopOnError) return kFALSE;
 505   }
 506
 507   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 508                stopwatch.RealTime(),stopwatch.CpuTime()));
 509
 510   return kTRUE;
 511 }
 512
 513 //_____________________________________________________________________________
 514 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 515 {
 516 // run the local reconstruction
 517
 518   TStopwatch stopwatch;
 519   stopwatch.Start();
 520
 521   TString detStr = detectors;
 522   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 523     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 524     AliReconstructor* reconstructor = GetReconstructor(iDet);
 525     if (!reconstructor) continue;
 526     AliLoader* loader = fLoader[iDet];
 527
 528     // conversion of digits
 529     if (fRawReader && reconstructor->HasDigitConversion()) {
 530       AliInfo(Form("converting raw data digits into root objects for %s",
 531                    fgkDetectorName[iDet]));
 532       TStopwatch stopwatchDet;
 533       stopwatchDet.Start();
 534       loader->LoadDigits("update");
 535       loader->CleanDigits();
 536       loader->MakeDigitsContainer();
 537       TTree* digitsTree = loader->TreeD();
 538       reconstructor->ConvertDigits(fRawReader, digitsTree);
 539       loader->WriteDigits("OVERWRITE");
 540       loader->UnloadDigits();
 541       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 542                    fgkDetectorName[iDet],
 543                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 544     }
 545
 546     // local reconstruction
 547     if (!reconstructor->HasLocalReconstruction()) continue;
 548     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 549     TStopwatch stopwatchDet;
 550     stopwatchDet.Start();
 551     loader->LoadRecPoints("update");
 552     loader->CleanRecPoints();
 553     loader->MakeRecPointsContainer();
 554     TTree* clustersTree = loader->TreeR();
 555     if (fRawReader && !reconstructor->HasDigitConversion()) {
 556       reconstructor->Reconstruct(fRawReader, clustersTree);
 557     } else {
 558       loader->LoadDigits("read");
 559       TTree* digitsTree = loader->TreeD();
 560       if (!digitsTree) {
 561         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 562         if (fStopOnError) return kFALSE;
 563       } else {
 564         reconstructor->Reconstruct(digitsTree, clustersTree);
 565       }
 566       loader->UnloadDigits();
 567     }
 568     loader->WriteRecPoints("OVERWRITE");
 569     loader->UnloadRecPoints();
 570     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
 571                     fgkDetectorName[iDet],
 572                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 573   }
 574
 575   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 576     AliError(Form("the following detectors were not found: %s",
 577                   detStr.Data()));
 578     if (fStopOnError) return kFALSE;
 579   }
 580
 581   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 582                stopwatch.RealTime(),stopwatch.CpuTime()));
 583
 584   return kTRUE;
 585 }
 586
 587 //_____________________________________________________________________________
 588 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 589 {
 590 // run the barrel tracking
 591
 592   TStopwatch stopwatch;
 593   stopwatch.Start();
 594
 595   AliESDVertex* vertex = NULL;
 596   Double_t vtxPos[3] = {0, 0, 0};
 597   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 598   TArrayF mcVertex(3);
 599   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 600     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 601     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 602   }
 603
 604   if (fVertexer) {
 605     AliInfo("running the ITS vertex finder");
 606     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 607     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 608     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 609     if(!vertex){
 610       AliWarning("Vertex not found");
 611       vertex = new AliESDVertex();
 612     }
 613     else {
 614       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 615     }
 616
 617   } else {
 618     AliInfo("getting the primary vertex from MC");
 619     vertex = new AliESDVertex(vtxPos, vtxErr);
 620   }
 621
 622   if (vertex) {
 623     vertex->GetXYZ(vtxPos);
 624     vertex->GetSigmaXYZ(vtxErr);
 625   } else {
 626     AliWarning("no vertex reconstructed");
 627     vertex = new AliESDVertex(vtxPos, vtxErr);
 628   }
 629   esd->SetVertex(vertex);
 630   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 631     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
 632   }
 633   delete vertex;
 634
 635   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 636                stopwatch.RealTime(),stopwatch.CpuTime()));
 637
 638   return kTRUE;
 639 }
 640
 641 //_____________________________________________________________________________
 642 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
 643 {
 644 // run the HLT barrel tracking
 645
 646   TStopwatch stopwatch;
 647   stopwatch.Start();
 648
 649   if (!fRunLoader) {
 650     AliError("Missing runLoader!");
 651     return kFALSE;
 652   }
 653
 654   AliInfo("running HLT tracking");
 655
 656   // Get a pointer to the HLT reconstructor
 657   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
 658   if (!reconstructor) return kFALSE;
 659
 660   // TPC + ITS
 661   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 662     TString detName = fgkDetectorName[iDet];
 663     AliDebug(1, Form("%s HLT tracking", detName.Data()));
 664     reconstructor->SetOption(detName.Data());
 665     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
 666     if (!tracker) {
 667       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
 668       if (fStopOnError) return kFALSE;
 669       continue;
 670     }
 671     Double_t vtxPos[3];
 672     Double_t vtxErr[3]={0.005,0.005,0.010};
 673     const AliESDVertex *vertex = esd->GetVertex();
 674     vertex->GetXYZ(vtxPos);
 675     tracker->SetVertex(vtxPos,vtxErr);
 676     if(iDet != 1) {
 677       fLoader[iDet]->LoadRecPoints("read");
 678       TTree* tree = fLoader[iDet]->TreeR();
 679       if (!tree) {
 680         AliError(Form("Can't get the %s cluster tree", detName.Data()));
 681         return kFALSE;
 682       }
 683       tracker->LoadClusters(tree);
 684     }
 685     if (tracker->Clusters2Tracks(esd) != 0) {
 686       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
 687       return kFALSE;
 688     }
 689     if(iDet != 1) {
 690       tracker->UnloadClusters();
 691     }
 692     delete tracker;
 693   }
 694
 695   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 696                stopwatch.RealTime(),stopwatch.CpuTime()));
 697
 698   return kTRUE;
 699 }
 700
 701 //_____________________________________________________________________________
 702 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
 703 {
 704 // run the barrel tracking
 705
 706   TStopwatch stopwatch;
 707   stopwatch.Start();
 708
 709   AliInfo("running tracking");
 710
 711   // pass 1: TPC + ITS inwards
 712   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 713     if (!fTracker[iDet]) continue;
 714     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
 715
 716     // load clusters
 717     fLoader[iDet]->LoadRecPoints("read");
 718     TTree* tree = fLoader[iDet]->TreeR();
 719     if (!tree) {
 720       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 721       return kFALSE;
 722     }
 723     fTracker[iDet]->LoadClusters(tree);
 724
 725     // run tracking
 726     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
 727       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
 728       return kFALSE;
 729     }
 730     if (fCheckPointLevel > 1) {
 731       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
 732     }
 733     // preliminary PID in TPC needed by the ITS tracker
 734     if (iDet == 1) {
 735       GetReconstructor(1)->FillESD(fRunLoader, esd);
 736       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
 737       AliESDpid::MakePID(esd);
 738     }
 739   }
 740
 741   // pass 2: ALL backwards
 742   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 743     if (!fTracker[iDet]) continue;
 744     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
 745
 746     // load clusters
 747     if (iDet > 1) {     // all except ITS, TPC
 748       TTree* tree = NULL;
 749       fLoader[iDet]->LoadRecPoints("read");
 750       tree = fLoader[iDet]->TreeR();
 751       if (!tree) {
 752         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 753         return kFALSE;
 754       }
 755       fTracker[iDet]->LoadClusters(tree);
 756     }
 757
 758     // run tracking
 759     if (fTracker[iDet]->PropagateBack(esd) != 0) {
 760       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
 761       return kFALSE;
 762     }
 763     if (fCheckPointLevel > 1) {
 764       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
 765     }
 766
 767     // unload clusters
 768     if (iDet > 2) {     // all except ITS, TPC, TRD
 769       fTracker[iDet]->UnloadClusters();
 770       fLoader[iDet]->UnloadRecPoints();
 771     }
 772     // updated PID in TPC needed by the ITS tracker -MI
 773     if (iDet == 1) {
 774       GetReconstructor(1)->FillESD(fRunLoader, esd);
 775       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
 776       AliESDpid::MakePID(esd);
 777     }
 778   }
 779
 780   // write space-points to the ESD in case alignment data output
 781   // is switched on
 782   if (fWriteAlignmentData)
 783     WriteAlignmentData(esd);
 784
 785   // pass 3: TRD + TPC + ITS refit inwards
 786   for (Int_t iDet = 2; iDet >= 0; iDet--) {
 787     if (!fTracker[iDet]) continue;
 788     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
 789
 790     // run tracking
 791     if (fTracker[iDet]->RefitInward(esd) != 0) {
 792       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
 793       return kFALSE;
 794     }
 795     if (fCheckPointLevel > 1) {
 796       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
 797     }
 798
 799     // unload clusters
 800     fTracker[iDet]->UnloadClusters();
 801     fLoader[iDet]->UnloadRecPoints();
 802   }
 803   //
 804   // Propagate track to the vertex - if not done by ITS
 805   //
 806   Int_t ntracks = esd->GetNumberOfTracks();
 807   for (Int_t itrack=0; itrack<ntracks; itrack++){
 808     const Double_t kRadius  = 3;   // beam pipe radius
 809     const Double_t kMaxStep = 5;   // max step
 810     const Double_t kMaxD    = 123456;  // max distance to prim vertex
 811     Double_t       fieldZ   = AliTracker::GetBz();  //
 812     AliESDtrack * track = esd->GetTrack(itrack);
 813     if (!track) continue;
 814     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
 815     track->PropagateTo(kRadius, track->GetMass(),kMaxStep,kTRUE);
 816     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
 817   }
 818
 819   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 820                stopwatch.RealTime(),stopwatch.CpuTime()));
 821
 822   return kTRUE;
 823 }
 824
 825 //_____________________________________________________________________________
 826 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
 827 {
 828 // fill the event summary data
 829
 830   TStopwatch stopwatch;
 831   stopwatch.Start();
 832   AliInfo("filling ESD");
 833
 834   TString detStr = detectors;
 835   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 836     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 837     AliReconstructor* reconstructor = GetReconstructor(iDet);
 838     if (!reconstructor) continue;
 839
 840     if (!ReadESD(esd, fgkDetectorName[iDet])) {
 841       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
 842       TTree* clustersTree = NULL;
 843       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 844         fLoader[iDet]->LoadRecPoints("read");
 845         clustersTree = fLoader[iDet]->TreeR();
 846         if (!clustersTree) {
 847           AliError(Form("Can't get the %s clusters tree",
 848                         fgkDetectorName[iDet]));
 849           if (fStopOnError) return kFALSE;
 850         }
 851       }
 852       if (fRawReader && !reconstructor->HasDigitConversion()) {
 853         reconstructor->FillESD(fRawReader, clustersTree, esd);
 854       } else {
 855         TTree* digitsTree = NULL;
 856         if (fLoader[iDet]) {
 857           fLoader[iDet]->LoadDigits("read");
 858           digitsTree = fLoader[iDet]->TreeD();
 859           if (!digitsTree) {
 860             AliError(Form("Can't get the %s digits tree",
 861                           fgkDetectorName[iDet]));
 862             if (fStopOnError) return kFALSE;
 863           }
 864         }
 865         reconstructor->FillESD(digitsTree, clustersTree, esd);
 866         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
 867       }
 868       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 869         fLoader[iDet]->UnloadRecPoints();
 870       }
 871
 872       if (fRawReader) {
 873         reconstructor->FillESD(fRunLoader, fRawReader, esd);
 874       } else {
 875         reconstructor->FillESD(fRunLoader, esd);
 876       }
 877       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
 878     }
 879   }
 880
 881   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 882     AliError(Form("the following detectors were not found: %s",
 883                   detStr.Data()));
 884     if (fStopOnError) return kFALSE;
 885   }
 886
 887   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 888                stopwatch.RealTime(),stopwatch.CpuTime()));
 889
 890   return kTRUE;
 891 }
 892
 893
 894 //_____________________________________________________________________________
 895 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
 896 {
 897 // check whether detName is contained in detectors
 898 // if yes, it is removed from detectors
 899
 900   // check if all detectors are selected
 901   if ((detectors.CompareTo("ALL") == 0) ||
 902       detectors.BeginsWith("ALL ") ||
 903       detectors.EndsWith(" ALL") ||
 904       detectors.Contains(" ALL ")) {
 905     detectors = "ALL";
 906     return kTRUE;
 907   }
 908
 909   // search for the given detector
 910   Bool_t result = kFALSE;
 911   if ((detectors.CompareTo(detName) == 0) ||
 912       detectors.BeginsWith(detName+" ") ||
 913       detectors.EndsWith(" "+detName) ||
 914       detectors.Contains(" "+detName+" ")) {
 915     detectors.ReplaceAll(detName, "");
 916     result = kTRUE;
 917   }
 918
 919   // clean up the detectors string
 920   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
 921   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
 922   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
 923
 924   return result;
 925 }
 926
 927 //_____________________________________________________________________________
 928 Bool_t AliReconstruction::InitRunLoader()
 929 {
 930 // get or create the run loader
 931
 932   if (gAlice) delete gAlice;
 933   gAlice = NULL;
 934
 935   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
 936     // load all base libraries to get the loader classes
 937     TString libs = gSystem->GetLibraries();
 938     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 939       TString detName = fgkDetectorName[iDet];
 940       if (detName == "HLT") continue;
 941       if (libs.Contains("lib" + detName + "base.so")) continue;
 942       gSystem->Load("lib" + detName + "base.so");
 943     }
 944     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
 945     if (!fRunLoader) {
 946       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
 947       CleanUp();
 948       return kFALSE;
 949     }
 950     fRunLoader->CdGAFile();
 951     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
 952       if (fRunLoader->LoadgAlice() == 0) {
 953         gAlice = fRunLoader->GetAliRun();
 954         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
 955       }
 956     }
 957     if (!gAlice && !fRawReader) {
 958       AliError(Form("no gAlice object found in file %s",
 959                     fGAliceFileName.Data()));
 960       CleanUp();
 961       return kFALSE;
 962     }
 963
 964   } else {               // galice.root does not exist
 965     if (!fRawReader) {
 966       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
 967       CleanUp();
 968       return kFALSE;
 969     }
 970     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
 971                                     AliConfig::GetDefaultEventFolderName(),
 972                                     "recreate");
 973     if (!fRunLoader) {
 974       AliError(Form("could not create run loader in file %s",
 975                     fGAliceFileName.Data()));
 976       CleanUp();
 977       return kFALSE;
 978     }
 979     fRunLoader->MakeTree("E");
 980     Int_t iEvent = 0;
 981     while (fRawReader->NextEvent()) {
 982       fRunLoader->SetEventNumber(iEvent);
 983       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
 984                                      iEvent, iEvent);
 985       fRunLoader->MakeTree("H");
 986       fRunLoader->TreeE()->Fill();
 987       iEvent++;
 988     }
 989     fRawReader->RewindEvents();
 990     fRunLoader->WriteHeader("OVERWRITE");
 991     fRunLoader->CdGAFile();
 992     fRunLoader->Write(0, TObject::kOverwrite);
 993 //    AliTracker::SetFieldMap(???);
 994   }
 995
 996   return kTRUE;
 997 }
 998
 999 //_____________________________________________________________________________
1000 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1001 {
1002 // get the reconstructor object and the loader for a detector
1003
1004   if (fReconstructor[iDet]) return fReconstructor[iDet];
1005
1006   // load the reconstructor object
1007   TPluginManager* pluginManager = gROOT->GetPluginManager();
1008   TString detName = fgkDetectorName[iDet];
1009   TString recName = "Ali" + detName + "Reconstructor";
1010   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1011
1012   if (detName == "HLT") {
1013     if (!gROOT->GetClass("AliLevel3")) {
1014       gSystem->Load("libAliL3Src.so");
1015       gSystem->Load("libAliL3Misc.so");
1016       gSystem->Load("libAliL3Hough.so");
1017       gSystem->Load("libAliL3Comp.so");
1018     }
1019   }
1020
1021   AliReconstructor* reconstructor = NULL;
1022   // first check if a plugin is defined for the reconstructor
1023   TPluginHandler* pluginHandler =
1024     pluginManager->FindHandler("AliReconstructor", detName);
1025   // if not, add a plugin for it
1026   if (!pluginHandler) {
1027     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1028     TString libs = gSystem->GetLibraries();
1029     if (libs.Contains("lib" + detName + "base.so") ||
1030         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1031       pluginManager->AddHandler("AliReconstructor", detName,
1032                                 recName, detName + "rec", recName + "()");
1033     } else {
1034       pluginManager->AddHandler("AliReconstructor", detName,
1035                                 recName, detName, recName + "()");
1036     }
1037     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1038   }
1039   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1040     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1041   }
1042   if (reconstructor) {
1043     TObject* obj = fOptions.FindObject(detName.Data());
1044     if (obj) reconstructor->SetOption(obj->GetTitle());
1045     reconstructor->Init(fRunLoader);
1046     fReconstructor[iDet] = reconstructor;
1047   }
1048
1049   // get or create the loader
1050   if (detName != "HLT") {
1051     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1052     if (!fLoader[iDet]) {
1053       AliConfig::Instance()
1054         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1055                                 detName, detName);
1056       // first check if a plugin is defined for the loader
1057       TPluginHandler* pluginHandler =
1058         pluginManager->FindHandler("AliLoader", detName);
1059       // if not, add a plugin for it
1060       if (!pluginHandler) {
1061         TString loaderName = "Ali" + detName + "Loader";
1062         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1063         pluginManager->AddHandler("AliLoader", detName,
1064                                   loaderName, detName + "base",
1065                                   loaderName + "(const char*, TFolder*)");
1066         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1067       }
1068       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1069         fLoader[iDet] =
1070           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1071                                                  fRunLoader->GetEventFolder());
1072       }
1073       if (!fLoader[iDet]) {   // use default loader
1074         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1075       }
1076       if (!fLoader[iDet]) {
1077         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1078         if (fStopOnError) return NULL;
1079       } else {
1080         fRunLoader->AddLoader(fLoader[iDet]);
1081         fRunLoader->CdGAFile();
1082         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1083         fRunLoader->Write(0, TObject::kOverwrite);
1084       }
1085     }
1086   }
1087
1088   return reconstructor;
1089 }
1090
1091 //_____________________________________________________________________________
1092 Bool_t AliReconstruction::CreateVertexer()
1093 {
1094 // create the vertexer
1095
1096   fVertexer = NULL;
1097   AliReconstructor* itsReconstructor = GetReconstructor(0);
1098   if (itsReconstructor) {
1099     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1100   }
1101   if (!fVertexer) {
1102     AliWarning("couldn't create a vertexer for ITS");
1103     if (fStopOnError) return kFALSE;
1104   }
1105
1106   return kTRUE;
1107 }
1108
1109 //_____________________________________________________________________________
1110 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1111 {
1112 // create the trackers
1113
1114   TString detStr = detectors;
1115   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1116     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1117     AliReconstructor* reconstructor = GetReconstructor(iDet);
1118     if (!reconstructor) continue;
1119     TString detName = fgkDetectorName[iDet];
1120     if (detName == "HLT") {
1121       fRunHLTTracking = kTRUE;
1122       continue;
1123     }
1124
1125     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1126     if (!fTracker[iDet] && (iDet < 7)) {
1127       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1128       if (fStopOnError) return kFALSE;
1129     }
1130   }
1131
1132   return kTRUE;
1133 }
1134
1135 //_____________________________________________________________________________
1136 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1137 {
1138 // delete trackers and the run loader and close and delete the file
1139
1140   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1141     delete fReconstructor[iDet];
1142     fReconstructor[iDet] = NULL;
1143     fLoader[iDet] = NULL;
1144     delete fTracker[iDet];
1145     fTracker[iDet] = NULL;
1146   }
1147   delete fVertexer;
1148   fVertexer = NULL;
1149
1150   delete fRunLoader;
1151   fRunLoader = NULL;
1152   delete fRawReader;
1153   fRawReader = NULL;
1154
1155   if (file) {
1156     file->Close();
1157     delete file;
1158   }
1159
1160   if (fileOld) {
1161     fileOld->Close();
1162     delete fileOld;
1163     gSystem->Unlink("AliESDs.old.root");
1164   }
1165 }
1166
1167
1168 //_____________________________________________________________________________
1169 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1170 {
1171 // read the ESD event from a file
1172
1173   if (!esd) return kFALSE;
1174   char fileName[256];
1175   sprintf(fileName, "ESD_%d.%d_%s.root",
1176           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1177   if (gSystem->AccessPathName(fileName)) return kFALSE;
1178
1179   AliInfo(Form("reading ESD from file %s", fileName));
1180   AliDebug(1, Form("reading ESD from file %s", fileName));
1181   TFile* file = TFile::Open(fileName);
1182   if (!file || !file->IsOpen()) {
1183     AliError(Form("opening %s failed", fileName));
1184     delete file;
1185     return kFALSE;
1186   }
1187
1188   gROOT->cd();
1189   delete esd;
1190   esd = (AliESD*) file->Get("ESD");
1191   file->Close();
1192   delete file;
1193   return kTRUE;
1194 }
1195
1196 //_____________________________________________________________________________
1197 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1198 {
1199 // write the ESD event to a file
1200
1201   if (!esd) return;
1202   char fileName[256];
1203   sprintf(fileName, "ESD_%d.%d_%s.root",
1204           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1205
1206   AliDebug(1, Form("writing ESD to file %s", fileName));
1207   TFile* file = TFile::Open(fileName, "recreate");
1208   if (!file || !file->IsOpen()) {
1209     AliError(Form("opening %s failed", fileName));
1210   } else {
1211     esd->Write("ESD");
1212     file->Close();
1213   }
1214   delete file;
1215 }
1216
1217
1218
1219
1220 //_____________________________________________________________________________
1221 void AliReconstruction::CreateTag(TFile* file)
1222 {
1223   /////////////
1224   //muon code//
1225   ////////////
1226   Double_t fMUONMASS = 0.105658369;
1227   //Variables
1228   Double_t fX,fY,fZ ;
1229   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1230   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1231   Int_t fCharge;
1232   TLorentzVector fEPvector;
1233
1234   Float_t fZVertexCut = 10.0;
1235   Float_t fRhoVertexCut = 2.0;
1236
1237   Float_t fLowPtCut = 1.0;
1238   Float_t fHighPtCut = 3.0;
1239   Float_t fVeryHighPtCut = 10.0;
1240   ////////////
1241
1242   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1243
1244   // Creates the tags for all the events in a given ESD file
1245   Int_t ntrack;
1246   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1247   Int_t nPos, nNeg, nNeutr;
1248   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1249   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1250   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1251   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1252   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1253
1254   AliRunTag *tag = new AliRunTag();
1255   AliEventTag *evTag = new AliEventTag();
1256   TTree ttag("T","A Tree with event tags");
1257   TBranch * btag = ttag.Branch("AliTAG", &tag);
1258   btag->SetCompressionLevel(9);
1259
1260   AliInfo(Form("Creating the tags......."));
1261
1262   if (!file || !file->IsOpen()) {
1263     AliError(Form("opening failed"));
1264     delete file;
1265     return ;
1266   }
1267   Int_t firstEvent = 0,lastEvent = 0;
1268   TTree *t = (TTree*) file->Get("esdTree");
1269   TBranch * b = t->GetBranch("ESD");
1270   AliESD *esd = 0;
1271   b->SetAddress(&esd);
1272
1273   tag->SetRunId(esd->GetRunNumber());
1274
1275   Int_t iNumberOfEvents = b->GetEntries();
1276   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1277     ntrack = 0;
1278     nPos = 0;
1279     nNeg = 0;
1280     nNeutr =0;
1281     nK0s = 0;
1282     nNeutrons = 0;
1283     nPi0s = 0;
1284     nGamas = 0;
1285     nProtons = 0;
1286     nKaons = 0;
1287     nPions = 0;
1288     nMuons = 0;
1289     nElectrons = 0;
1290     nCh1GeV = 0;
1291     nCh3GeV = 0;
1292     nCh10GeV = 0;
1293     nMu1GeV = 0;
1294     nMu3GeV = 0;
1295     nMu10GeV = 0;
1296     nEl1GeV = 0;
1297     nEl3GeV = 0;
1298     nEl10GeV = 0;
1299     maxPt = .0;
1300     meanPt = .0;
1301     totalP = .0;
1302
1303     b->GetEntry(iEventNumber);
1304     const AliESDVertex * vertexIn = esd->GetVertex();
1305
1306     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1307       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1308       UInt_t status = esdTrack->GetStatus();
1309
1310       //select only tracks with ITS refit
1311       if ((status&AliESDtrack::kITSrefit)==0) continue;
1312       //select only tracks with TPC refit
1313       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1314
1315       //select only tracks with the "combined PID"
1316       if ((status&AliESDtrack::kESDpid)==0) continue;
1317       Double_t p[3];
1318       esdTrack->GetPxPyPz(p);
1319       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1320       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1321       totalP += momentum;
1322       meanPt += fPt;
1323       if(fPt > maxPt) maxPt = fPt;
1324
1325       if(esdTrack->GetSign() > 0) {
1326         nPos++;
1327         if(fPt > fLowPtCut) nCh1GeV++;
1328         if(fPt > fHighPtCut) nCh3GeV++;
1329         if(fPt > fVeryHighPtCut) nCh10GeV++;
1330       }
1331       if(esdTrack->GetSign() < 0) {
1332         nNeg++;
1333         if(fPt > fLowPtCut) nCh1GeV++;
1334         if(fPt > fHighPtCut) nCh3GeV++;
1335         if(fPt > fVeryHighPtCut) nCh10GeV++;
1336       }
1337       if(esdTrack->GetSign() == 0) nNeutr++;
1338
1339       //PID
1340       Double_t prob[5];
1341       esdTrack->GetESDpid(prob);
1342
1343       Double_t rcc = 0.0;
1344       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1345       if(rcc == 0.0) continue;
1346       //Bayes' formula
1347       Double_t w[5];
1348       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1349
1350       //protons
1351       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1352       //kaons
1353       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1354       //pions
1355       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1356       //electrons
1357       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1358         nElectrons++;
1359         if(fPt > fLowPtCut) nEl1GeV++;
1360         if(fPt > fHighPtCut) nEl3GeV++;
1361         if(fPt > fVeryHighPtCut) nEl10GeV++;
1362       }
1363       ntrack++;
1364     }//track loop
1365
1366     /////////////
1367     //muon code//
1368     ////////////
1369     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1370     // loop over all reconstructed tracks (also first track of combination)
1371     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1372       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1373       if (muonTrack == 0x0) continue;
1374
1375       // Coordinates at vertex
1376       fZ = muonTrack->GetZ();
1377       fY = muonTrack->GetBendingCoor();
1378       fX = muonTrack->GetNonBendingCoor();
1379
1380       fThetaX = muonTrack->GetThetaX();
1381       fThetaY = muonTrack->GetThetaY();
1382
1383       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1384       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1385       fPxRec = fPzRec * TMath::Tan(fThetaX);
1386       fPyRec = fPzRec * TMath::Tan(fThetaY);
1387       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1388
1389       //ChiSquare of the track if needed
1390       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1391       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1392       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1393
1394       // total number of muons inside a vertex cut
1395       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1396         nMuons++;
1397         if(fEPvector.Pt() > fLowPtCut) {
1398           nMu1GeV++;
1399           if(fEPvector.Pt() > fHighPtCut) {
1400             nMu3GeV++;
1401             if (fEPvector.Pt() > fVeryHighPtCut) {
1402               nMu10GeV++;
1403             }
1404           }
1405         }
1406       }
1407     }//muon track loop
1408
1409     // Fill the event tags
1410     if(ntrack != 0)
1411       meanPt = meanPt/ntrack;
1412
1413     evTag->SetEventId(iEventNumber+1);
1414     evTag->SetVertexX(vertexIn->GetXv());
1415     evTag->SetVertexY(vertexIn->GetYv());
1416     evTag->SetVertexZ(vertexIn->GetZv());
1417
1418     evTag->SetT0VertexZ(esd->GetT0zVertex());
1419
1420     evTag->SetTrigger(esd->GetTrigger());
1421
1422     evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
1423     evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
1424     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1425     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1426
1427
1428     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1429     evTag->SetNumOfPosTracks(nPos);
1430     evTag->SetNumOfNegTracks(nNeg);
1431     evTag->SetNumOfNeutrTracks(nNeutr);
1432
1433     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1434     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1435     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1436     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1437
1438     evTag->SetNumOfProtons(nProtons);
1439     evTag->SetNumOfKaons(nKaons);
1440     evTag->SetNumOfPions(nPions);
1441     evTag->SetNumOfMuons(nMuons);
1442     evTag->SetNumOfElectrons(nElectrons);
1443     evTag->SetNumOfPhotons(nGamas);
1444     evTag->SetNumOfPi0s(nPi0s);
1445     evTag->SetNumOfNeutrons(nNeutrons);
1446     evTag->SetNumOfKaon0s(nK0s);
1447
1448     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1449     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1450     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1451     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1452     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1453     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1454     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1455     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1456     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1457
1458     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1459     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1460
1461     evTag->SetTotalMomentum(totalP);
1462     evTag->SetMeanPt(meanPt);
1463     evTag->SetMaxPt(maxPt);
1464
1465     tag->AddEventTag(*evTag);
1466   }
1467   lastEvent = iNumberOfEvents;
1468
1469   ttag.Fill();
1470   tag->Clear();
1471
1472   char fileName[256];
1473   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1474           tag->GetRunId(),firstEvent,lastEvent );
1475   AliInfo(Form("writing tags to file %s", fileName));
1476   AliDebug(1, Form("writing tags to file %s", fileName));
1477
1478   TFile* ftag = TFile::Open(fileName, "recreate");
1479   ftag->cd();
1480   ttag.Write();
1481   ftag->Close();
1482   file->cd();
1483   delete tag;
1484   delete evTag;
1485 }
1486
1487 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1488 {
1489   // Write space-points which are then used in the alignment procedures
1490   // For the moment only ITS, TRD and TPC
1491
1492   // Load TOF clusters
1493   if (fTracker[3]){
1494     fLoader[3]->LoadRecPoints("read");
1495     TTree* tree = fLoader[3]->TreeR();
1496     if (!tree) {
1497       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1498       return;
1499     }
1500     fTracker[3]->LoadClusters(tree);
1501   }
1502   Int_t ntracks = esd->GetNumberOfTracks();
1503   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1504     {
1505       AliESDtrack *track = esd->GetTrack(itrack);
1506       Int_t nsp = 0;
1507       UInt_t idx[200];
1508       for (Int_t iDet = 3; iDet >= 0; iDet--)
1509         nsp += track->GetNcls(iDet);
1510       if (nsp) {
1511         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1512         track->SetTrackPointArray(sp);
1513         Int_t isptrack = 0;
1514         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1515           AliTracker *tracker = fTracker[iDet];
1516           if (!tracker) continue;
1517           Int_t nspdet = track->GetNcls(iDet);
1518           if (nspdet <= 0) continue;
1519           track->GetClusters(iDet,idx);
1520           AliTrackPoint p;
1521           Int_t isp = 0;
1522           Int_t isp2 = 0;
1523           while (isp < nspdet) {
1524             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1525             if (!isvalid) continue;
1526             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1527           }
1528         }
1529       }
1530     }
1531   if (fTracker[3]){
1532     fTracker[3]->UnloadClusters();
1533     fLoader[3]->UnloadRecPoints();
1534   }
1535 }