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