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