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