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       }
 911     }
 912     if (!gAlice && !fRawReader) {
 913       AliError(Form("no gAlice object found in file %s",
 914                     fGAliceFileName.Data()));
 915       CleanUp();
 916       return kFALSE;
 917     }
 918
 919   } else {               // galice.root does not exist
 920     if (!fRawReader) {
 921       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
 922       CleanUp();
 923       return kFALSE;
 924     }
 925     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
 926                                     AliConfig::GetDefaultEventFolderName(),
 927                                     "recreate");
 928     if (!fRunLoader) {
 929       AliError(Form("could not create run loader in file %s",
 930                     fGAliceFileName.Data()));
 931       CleanUp();
 932       return kFALSE;
 933     }
 934     fRunLoader->MakeTree("E");
 935     Int_t iEvent = 0;
 936     while (fRawReader->NextEvent()) {
 937       fRunLoader->SetEventNumber(iEvent);
 938       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
 939                                      iEvent, iEvent);
 940       fRunLoader->MakeTree("H");
 941       fRunLoader->TreeE()->Fill();
 942       iEvent++;
 943     }
 944     fRawReader->RewindEvents();
 945     fRunLoader->WriteHeader("OVERWRITE");
 946     fRunLoader->CdGAFile();
 947     fRunLoader->Write(0, TObject::kOverwrite);
 948 //    AliTracker::SetFieldMap(???);
 949   }
 950
 951   return kTRUE;
 952 }
 953
 954 //_____________________________________________________________________________
 955 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
 956 {
 957 // get the reconstructor object and the loader for a detector
 958
 959   if (fReconstructor[iDet]) return fReconstructor[iDet];
 960
 961   // load the reconstructor object
 962   TPluginManager* pluginManager = gROOT->GetPluginManager();
 963   TString detName = fgkDetectorName[iDet];
 964   TString recName = "Ali" + detName + "Reconstructor";
 965   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
 966
 967   if (detName == "HLT") {
 968     if (!gROOT->GetClass("AliLevel3")) {
 969       gSystem->Load("libAliL3Src.so");
 970       gSystem->Load("libAliL3Misc.so");
 971       gSystem->Load("libAliL3Hough.so");
 972       gSystem->Load("libAliL3Comp.so");
 973     }
 974   }
 975
 976   AliReconstructor* reconstructor = NULL;
 977   // first check if a plugin is defined for the reconstructor
 978   TPluginHandler* pluginHandler =
 979     pluginManager->FindHandler("AliReconstructor", detName);
 980   // if not, add a plugin for it
 981   if (!pluginHandler) {
 982     AliDebug(1, Form("defining plugin for %s", recName.Data()));
 983     TString libs = gSystem->GetLibraries();
 984     if (libs.Contains("lib" + detName + "base.so") ||
 985         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
 986       pluginManager->AddHandler("AliReconstructor", detName,
 987                                 recName, detName + "rec", recName + "()");
 988     } else {
 989       pluginManager->AddHandler("AliReconstructor", detName,
 990                                 recName, detName, recName + "()");
 991     }
 992     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
 993   }
 994   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
 995     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
 996   }
 997   if (reconstructor) {
 998     TObject* obj = fOptions.FindObject(detName.Data());
 999     if (obj) reconstructor->SetOption(obj->GetTitle());
1000     reconstructor->Init(fRunLoader);
1001     fReconstructor[iDet] = reconstructor;
1002   }
1003
1004   // get or create the loader
1005   if (detName != "HLT") {
1006     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1007     if (!fLoader[iDet]) {
1008       AliConfig::Instance()
1009         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1010                                 detName, detName);
1011       // first check if a plugin is defined for the loader
1012       TPluginHandler* pluginHandler =
1013         pluginManager->FindHandler("AliLoader", detName);
1014       // if not, add a plugin for it
1015       if (!pluginHandler) {
1016         TString loaderName = "Ali" + detName + "Loader";
1017         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1018         pluginManager->AddHandler("AliLoader", detName,
1019                                   loaderName, detName + "base",
1020                                   loaderName + "(const char*, TFolder*)");
1021         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1022       }
1023       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1024         fLoader[iDet] =
1025           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1026                                                  fRunLoader->GetEventFolder());
1027       }
1028       if (!fLoader[iDet]) {   // use default loader
1029         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1030       }
1031       if (!fLoader[iDet]) {
1032         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1033         if (fStopOnError) return NULL;
1034       } else {
1035         fRunLoader->AddLoader(fLoader[iDet]);
1036         fRunLoader->CdGAFile();
1037         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1038         fRunLoader->Write(0, TObject::kOverwrite);
1039       }
1040     }
1041   }
1042
1043   return reconstructor;
1044 }
1045
1046 //_____________________________________________________________________________
1047 Bool_t AliReconstruction::CreateVertexer()
1048 {
1049 // create the vertexer
1050
1051   fVertexer = NULL;
1052   AliReconstructor* itsReconstructor = GetReconstructor(0);
1053   if (itsReconstructor) {
1054     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1055   }
1056   if (!fVertexer) {
1057     AliWarning("couldn't create a vertexer for ITS");
1058     if (fStopOnError) return kFALSE;
1059   }
1060
1061   return kTRUE;
1062 }
1063
1064 //_____________________________________________________________________________
1065 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1066 {
1067 // create the trackers
1068
1069   TString detStr = detectors;
1070   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1071     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1072     AliReconstructor* reconstructor = GetReconstructor(iDet);
1073     if (!reconstructor) continue;
1074     TString detName = fgkDetectorName[iDet];
1075     if (detName == "HLT") {
1076       fRunHLTTracking = kTRUE;
1077       continue;
1078     }
1079
1080     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1081     if (!fTracker[iDet] && (iDet < 7)) {
1082       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1083       if (fStopOnError) return kFALSE;
1084     }
1085   }
1086
1087   return kTRUE;
1088 }
1089
1090 //_____________________________________________________________________________
1091 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1092 {
1093 // delete trackers and the run loader and close and delete the file
1094
1095   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1096     delete fReconstructor[iDet];
1097     fReconstructor[iDet] = NULL;
1098     fLoader[iDet] = NULL;
1099     delete fTracker[iDet];
1100     fTracker[iDet] = NULL;
1101   }
1102   delete fVertexer;
1103   fVertexer = NULL;
1104
1105   delete fRunLoader;
1106   fRunLoader = NULL;
1107   delete fRawReader;
1108   fRawReader = NULL;
1109
1110   if (file) {
1111     file->Close();
1112     delete file;
1113   }
1114
1115   if (fileOld) {
1116     fileOld->Close();
1117     delete fileOld;
1118     gSystem->Unlink("AliESDs.old.root");
1119   }
1120 }
1121
1122
1123 //_____________________________________________________________________________
1124 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1125 {
1126 // read the ESD event from a file
1127
1128   if (!esd) return kFALSE;
1129   char fileName[256];
1130   sprintf(fileName, "ESD_%d.%d_%s.root",
1131           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1132   if (gSystem->AccessPathName(fileName)) return kFALSE;
1133
1134   AliInfo(Form("reading ESD from file %s", fileName));
1135   AliDebug(1, Form("reading ESD from file %s", fileName));
1136   TFile* file = TFile::Open(fileName);
1137   if (!file || !file->IsOpen()) {
1138     AliError(Form("opening %s failed", fileName));
1139     delete file;
1140     return kFALSE;
1141   }
1142
1143   gROOT->cd();
1144   delete esd;
1145   esd = (AliESD*) file->Get("ESD");
1146   file->Close();
1147   delete file;
1148   return kTRUE;
1149 }
1150
1151 //_____________________________________________________________________________
1152 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1153 {
1154 // write the ESD event to a file
1155
1156   if (!esd) return;
1157   char fileName[256];
1158   sprintf(fileName, "ESD_%d.%d_%s.root",
1159           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1160
1161   AliDebug(1, Form("writing ESD to file %s", fileName));
1162   TFile* file = TFile::Open(fileName, "recreate");
1163   if (!file || !file->IsOpen()) {
1164     AliError(Form("opening %s failed", fileName));
1165   } else {
1166     esd->Write("ESD");
1167     file->Close();
1168   }
1169   delete file;
1170 }
1171
1172
1173
1174
1175 //_____________________________________________________________________________
1176 void AliReconstruction::CreateTag(TFile* file)
1177 {
1178   Int_t ntrack;
1179   Int_t NProtons, NKaons, NPions, NMuons, NElectrons;
1180   Int_t Npos, Nneg, Nneutr;
1181   Int_t NK0s, Nneutrons, Npi0s, Ngamas;
1182   Int_t Nch1GeV, Nch3GeV, Nch10GeV;
1183   Int_t Nmu1GeV, Nmu3GeV, Nmu10GeV;
1184   Int_t Nel1GeV, Nel3GeV, Nel10GeV;
1185   Float_t MaxPt = .0, MeanPt = .0, TotalP = .0;
1186
1187   AliRunTag *tag = new AliRunTag();
1188   AliDetectorTag *detTag = new AliDetectorTag();
1189   AliEventTag *evTag = new AliEventTag();
1190   TTree ttag("T","A Tree with event tags");
1191   TBranch * btag = ttag.Branch("AliTAG", "AliRunTag", &tag);
1192   btag->SetCompressionLevel(9);
1193
1194   AliInfo(Form("Creating the tags......."));
1195
1196   if (!file || !file->IsOpen()) {
1197     AliError(Form("opening failed"));
1198     delete file;
1199     return ;
1200   }
1201
1202   TTree *t = (TTree*) file->Get("esdTree");
1203   TBranch * b = t->GetBranch("ESD");
1204   AliESD *esd = 0;
1205   b->SetAddress(&esd);
1206
1207   tag->SetRunId(esd->GetRunNumber());
1208
1209   Int_t firstEvent = 0,lastEvent = 0;
1210   Int_t i_NumberOfEvents = b->GetEntries();
1211   for (Int_t i_EventNumber = 0; i_EventNumber < i_NumberOfEvents; i_EventNumber++)
1212     {
1213       ntrack = 0;
1214       Npos = 0;
1215       Nneg = 0;
1216       Nneutr =0;
1217       NK0s = 0;
1218       Nneutrons = 0;
1219       Npi0s = 0;
1220       Ngamas = 0;
1221       NProtons = 0;
1222       NKaons = 0;
1223       NPions = 0;
1224       NMuons = 0;
1225       NElectrons = 0;
1226       Nch1GeV = 0;
1227       Nch3GeV = 0;
1228       Nch10GeV = 0;
1229       Nmu1GeV = 0;
1230       Nmu3GeV = 0;
1231       Nmu10GeV = 0;
1232       Nel1GeV = 0;
1233       Nel3GeV = 0;
1234       Nel10GeV = 0;
1235       MaxPt = .0;
1236       MeanPt = .0;
1237       TotalP = .0;
1238
1239       b->GetEntry(i_EventNumber);
1240       const AliESDVertex * VertexIn = esd->GetVertex();
1241
1242       for (Int_t i_TrackNumber = 0; i_TrackNumber < esd->GetNumberOfTracks(); i_TrackNumber++)
1243         {
1244           AliESDtrack * ESDTrack = esd->GetTrack(i_TrackNumber);
1245           UInt_t status = ESDTrack->GetStatus();
1246
1247           //select only tracks with ITS refit
1248           if ((status&AliESDtrack::kITSrefit)==0) continue;
1249
1250           //select only tracks with TPC refit-->remove extremely high Pt tracks
1251           if ((status&AliESDtrack::kTPCrefit)==0) continue;
1252
1253           //select only tracks with the "combined PID"
1254           if ((status&AliESDtrack::kESDpid)==0) continue;
1255                   Double_t p[3];
1256           ESDTrack->GetPxPyPz(p);
1257           Double_t P = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1258           Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1259           TotalP += P;
1260           MeanPt += fPt;
1261           if(fPt > MaxPt)
1262             MaxPt = fPt;
1263
1264           if(ESDTrack->GetSign() > 0)
1265             {
1266               Npos++;
1267               if(fPt > 1.0)
1268                 Nch1GeV++;
1269               if(fPt > 3.0)
1270                 Nch3GeV++;
1271               if(fPt > 10.0)
1272                 Nch10GeV++;
1273             }
1274           if(ESDTrack->GetSign() < 0)
1275             {
1276               Nneg++;
1277               if(fPt > 1.0)
1278                 Nch1GeV++;
1279               if(fPt > 3.0)
1280                 Nch3GeV++;
1281               if(fPt > 10.0)
1282                 Nch10GeV++;
1283             }
1284           if(ESDTrack->GetSign() == 0)
1285             Nneutr++;
1286
1287           //PID
1288           Double_t prob[10];
1289           ESDTrack->GetESDpid(prob);
1290
1291           //K0s
1292           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]))
1293             NK0s++;
1294           //neutrons
1295           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]))
1296             Nneutrons++;
1297           //pi0s
1298           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]))
1299             Npi0s++;
1300           //gamas
1301           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]))
1302             Ngamas++;
1303           //protons
1304           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]))
1305             NProtons++;
1306           //kaons
1307           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]))
1308             NKaons++;
1309           //kaons
1310           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]))
1311             NPions++;
1312           //muons
1313           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]))
1314             {
1315               NMuons++;
1316               if(fPt > 1.0)
1317                 Nmu1GeV++;
1318               if(fPt > 3.0)
1319                 Nmu3GeV++;
1320               if(fPt > 10.0)
1321                 Nmu10GeV++;
1322             }
1323           //electrons
1324           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]))
1325             {
1326               NElectrons++;
1327               if(fPt > 1.0)
1328                 Nel1GeV++;
1329               if(fPt > 3.0)
1330                 Nel3GeV++;
1331               if(fPt > 10.0)
1332                 Nel10GeV++;
1333             }
1334
1335
1336
1337           ntrack++;
1338         }//track loop
1339       // Fill the event tags
1340       MeanPt = MeanPt/ntrack;
1341
1342       evTag->SetEventId(i_EventNumber+1);
1343       evTag->SetVertexX(VertexIn->GetXv());
1344       evTag->SetVertexY(VertexIn->GetYv());
1345       evTag->SetVertexZ(VertexIn->GetZv());
1346
1347       evTag->SetT0VertexZ(esd->GetT0zVertex());
1348
1349       evTag->SetTrigger(esd->GetTrigger());
1350
1351       evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
1352       evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
1353       evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1354       evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1355
1356
1357       evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1358       evTag->SetNumOfPosTracks(Npos);
1359       evTag->SetNumOfNegTracks(Nneg);
1360       evTag->SetNumOfNeutrTracks(Nneutr);
1361
1362       evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1363       evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1364       evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1365       evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1366
1367       evTag->SetNumOfProtons(NProtons);
1368       evTag->SetNumOfKaons(NKaons);
1369       evTag->SetNumOfPions(NPions);
1370       evTag->SetNumOfMuons(NMuons);
1371       evTag->SetNumOfElectrons(NElectrons);
1372       evTag->SetNumOfPhotons(Ngamas);
1373       evTag->SetNumOfPi0s(Npi0s);
1374       evTag->SetNumOfNeutrons(Nneutrons);
1375       evTag->SetNumOfKaon0s(NK0s);
1376
1377       evTag->SetNumOfChargedAbove1GeV(Nch1GeV);
1378       evTag->SetNumOfChargedAbove3GeV(Nch3GeV);
1379       evTag->SetNumOfChargedAbove10GeV(Nch10GeV);
1380       evTag->SetNumOfMuonsAbove1GeV(Nmu1GeV);
1381       evTag->SetNumOfMuonsAbove3GeV(Nmu3GeV);
1382       evTag->SetNumOfMuonsAbove10GeV(Nmu10GeV);
1383       evTag->SetNumOfElectronsAbove1GeV(Nel1GeV);
1384       evTag->SetNumOfElectronsAbove3GeV(Nel3GeV);
1385       evTag->SetNumOfElectronsAbove10GeV(Nel10GeV);
1386
1387       evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
1388       evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
1389
1390       evTag->SetTotalMomentum(TotalP);
1391       evTag->SetMeanPt(MeanPt);
1392       evTag->SetMaxPt(MaxPt);
1393
1394       tag->AddEventTag(evTag);
1395     }
1396   lastEvent = i_NumberOfEvents;
1397
1398   ttag.Fill();
1399   tag->Clear();
1400
1401   char fileName[256];
1402   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1403           tag->GetRunId(),firstEvent,lastEvent );
1404   AliInfo(Form("writing tags to file %s", fileName));
1405   AliDebug(1, Form("writing tags to file %s", fileName));
1406
1407   TFile* ftag = TFile::Open(fileName, "recreate");
1408   ftag->cd();
1409   ttag.Write();
1410   ftag->Close();
1411   file->cd();
1412   delete tag;
1413   delete detTag;
1414   delete evTag;
1415 }
1416