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       fLoader[iDet]->LoadRecPoints("read");
 722       tree = fLoader[iDet]->TreeR();
 723       if (!tree) {
 724         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 725         return kFALSE;
 726       }
 727       fTracker[iDet]->LoadClusters(tree);
 728     }
 729
 730     // run tracking
 731     if (fTracker[iDet]->PropagateBack(esd) != 0) {
 732       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
 733       return kFALSE;
 734     }
 735     if (fCheckPointLevel > 1) {
 736       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
 737     }
 738
 739     // unload clusters
 740     if (iDet > 2) {     // all except ITS, TPC, TRD
 741       fTracker[iDet]->UnloadClusters();
 742       fLoader[iDet]->UnloadRecPoints();
 743     }
 744   }
 745
 746   // pass 3: TRD + TPC + ITS refit inwards
 747   for (Int_t iDet = 2; iDet >= 0; iDet--) {
 748     if (!fTracker[iDet]) continue;
 749     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
 750
 751     // run tracking
 752     if (fTracker[iDet]->RefitInward(esd) != 0) {
 753       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
 754       return kFALSE;
 755     }
 756     if (fCheckPointLevel > 1) {
 757       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
 758     }
 759
 760     // unload clusters
 761     fTracker[iDet]->UnloadClusters();
 762     fLoader[iDet]->UnloadRecPoints();
 763   }
 764
 765   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 766                stopwatch.RealTime(),stopwatch.CpuTime()));
 767
 768   return kTRUE;
 769 }
 770
 771 //_____________________________________________________________________________
 772 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
 773 {
 774 // fill the event summary data
 775
 776   TStopwatch stopwatch;
 777   stopwatch.Start();
 778   AliInfo("filling ESD");
 779
 780   TString detStr = detectors;
 781   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 782     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 783     AliReconstructor* reconstructor = GetReconstructor(iDet);
 784     if (!reconstructor) continue;
 785
 786     if (!ReadESD(esd, fgkDetectorName[iDet])) {
 787       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
 788       TTree* clustersTree = NULL;
 789       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 790         fLoader[iDet]->LoadRecPoints("read");
 791         clustersTree = fLoader[iDet]->TreeR();
 792         if (!clustersTree) {
 793           AliError(Form("Can't get the %s clusters tree",
 794                         fgkDetectorName[iDet]));
 795           if (fStopOnError) return kFALSE;
 796         }
 797       }
 798       if (fRawReader && !reconstructor->HasDigitConversion()) {
 799         reconstructor->FillESD(fRawReader, clustersTree, esd);
 800       } else {
 801         TTree* digitsTree = NULL;
 802         if (fLoader[iDet]) {
 803           fLoader[iDet]->LoadDigits("read");
 804           digitsTree = fLoader[iDet]->TreeD();
 805           if (!digitsTree) {
 806             AliError(Form("Can't get the %s digits tree",
 807                           fgkDetectorName[iDet]));
 808             if (fStopOnError) return kFALSE;
 809           }
 810         }
 811         reconstructor->FillESD(digitsTree, clustersTree, esd);
 812         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
 813       }
 814       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 815         fLoader[iDet]->UnloadRecPoints();
 816       }
 817
 818       if (fRawReader) {
 819         reconstructor->FillESD(fRunLoader, fRawReader, esd);
 820       } else {
 821         reconstructor->FillESD(fRunLoader, esd);
 822       }
 823       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
 824     }
 825   }
 826
 827   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 828     AliError(Form("the following detectors were not found: %s",
 829                   detStr.Data()));
 830     if (fStopOnError) return kFALSE;
 831   }
 832
 833   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 834                stopwatch.RealTime(),stopwatch.CpuTime()));
 835
 836   return kTRUE;
 837 }
 838
 839
 840 //_____________________________________________________________________________
 841 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
 842 {
 843 // check whether detName is contained in detectors
 844 // if yes, it is removed from detectors
 845
 846   // check if all detectors are selected
 847   if ((detectors.CompareTo("ALL") == 0) ||
 848       detectors.BeginsWith("ALL ") ||
 849       detectors.EndsWith(" ALL") ||
 850       detectors.Contains(" ALL ")) {
 851     detectors = "ALL";
 852     return kTRUE;
 853   }
 854
 855   // search for the given detector
 856   Bool_t result = kFALSE;
 857   if ((detectors.CompareTo(detName) == 0) ||
 858       detectors.BeginsWith(detName+" ") ||
 859       detectors.EndsWith(" "+detName) ||
 860       detectors.Contains(" "+detName+" ")) {
 861     detectors.ReplaceAll(detName, "");
 862     result = kTRUE;
 863   }
 864
 865   // clean up the detectors string
 866   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
 867   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
 868   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
 869
 870   return result;
 871 }
 872
 873 //_____________________________________________________________________________
 874 Bool_t AliReconstruction::InitRunLoader()
 875 {
 876 // get or create the run loader
 877
 878   if (gAlice) delete gAlice;
 879   gAlice = NULL;
 880
 881   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
 882     // load all base libraries to get the loader classes
 883     TString libs = gSystem->GetLibraries();
 884     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 885       TString detName = fgkDetectorName[iDet];
 886       if (detName == "HLT") continue;
 887       if (libs.Contains("lib" + detName + "base.so")) continue;
 888       gSystem->Load("lib" + detName + "base.so");
 889     }
 890     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
 891     if (!fRunLoader) {
 892       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
 893       CleanUp();
 894       return kFALSE;
 895     }
 896     fRunLoader->CdGAFile();
 897     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
 898       if (fRunLoader->LoadgAlice() == 0) {
 899         gAlice = fRunLoader->GetAliRun();
 900         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
 901         AliExternalTrackParam::SetFieldMap(gAlice->Field());
 902         if(fUniformField)
 903           AliExternalTrackParam::SetUniformFieldTracking();
 904         else
 905           AliExternalTrackParam::SetNonuniformFieldTracking();
 906       }
 907     }
 908     if (!gAlice && !fRawReader) {
 909       AliError(Form("no gAlice object found in file %s",
 910                     fGAliceFileName.Data()));
 911       CleanUp();
 912       return kFALSE;
 913     }
 914
 915   } else {               // galice.root does not exist
 916     if (!fRawReader) {
 917       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
 918       CleanUp();
 919       return kFALSE;
 920     }
 921     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
 922                                     AliConfig::GetDefaultEventFolderName(),
 923                                     "recreate");
 924     if (!fRunLoader) {
 925       AliError(Form("could not create run loader in file %s",
 926                     fGAliceFileName.Data()));
 927       CleanUp();
 928       return kFALSE;
 929     }
 930     fRunLoader->MakeTree("E");
 931     Int_t iEvent = 0;
 932     while (fRawReader->NextEvent()) {
 933       fRunLoader->SetEventNumber(iEvent);
 934       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
 935                                      iEvent, iEvent);
 936       fRunLoader->MakeTree("H");
 937       fRunLoader->TreeE()->Fill();
 938       iEvent++;
 939     }
 940     fRawReader->RewindEvents();
 941     fRunLoader->WriteHeader("OVERWRITE");
 942     fRunLoader->CdGAFile();
 943     fRunLoader->Write(0, TObject::kOverwrite);
 944 //    AliTracker::SetFieldMap(???);
 945   }
 946
 947   return kTRUE;
 948 }
 949
 950 //_____________________________________________________________________________
 951 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
 952 {
 953 // get the reconstructor object and the loader for a detector
 954
 955   if (fReconstructor[iDet]) return fReconstructor[iDet];
 956
 957   // load the reconstructor object
 958   TPluginManager* pluginManager = gROOT->GetPluginManager();
 959   TString detName = fgkDetectorName[iDet];
 960   TString recName = "Ali" + detName + "Reconstructor";
 961   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
 962
 963   if (detName == "HLT") {
 964     if (!gROOT->GetClass("AliLevel3")) {
 965       gSystem->Load("libAliL3Src.so");
 966       gSystem->Load("libAliL3Misc.so");
 967       gSystem->Load("libAliL3Hough.so");
 968       gSystem->Load("libAliL3Comp.so");
 969     }
 970   }
 971
 972   AliReconstructor* reconstructor = NULL;
 973   // first check if a plugin is defined for the reconstructor
 974   TPluginHandler* pluginHandler =
 975     pluginManager->FindHandler("AliReconstructor", detName);
 976   // if not, add a plugin for it
 977   if (!pluginHandler) {
 978     AliDebug(1, Form("defining plugin for %s", recName.Data()));
 979     TString libs = gSystem->GetLibraries();
 980     if (libs.Contains("lib" + detName + "base.so") ||
 981         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
 982       pluginManager->AddHandler("AliReconstructor", detName,
 983                                 recName, detName + "rec", recName + "()");
 984     } else {
 985       pluginManager->AddHandler("AliReconstructor", detName,
 986                                 recName, detName, recName + "()");
 987     }
 988     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
 989   }
 990   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
 991     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
 992   }
 993   if (reconstructor) {
 994     TObject* obj = fOptions.FindObject(detName.Data());
 995     if (obj) reconstructor->SetOption(obj->GetTitle());
 996     reconstructor->Init(fRunLoader);
 997     fReconstructor[iDet] = reconstructor;
 998   }
 999
1000   // get or create the loader
1001   if (detName != "HLT") {
1002     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1003     if (!fLoader[iDet]) {
1004       AliConfig::Instance()
1005         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1006                                 detName, detName);
1007       // first check if a plugin is defined for the loader
1008       TPluginHandler* pluginHandler =
1009         pluginManager->FindHandler("AliLoader", detName);
1010       // if not, add a plugin for it
1011       if (!pluginHandler) {
1012         TString loaderName = "Ali" + detName + "Loader";
1013         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1014         pluginManager->AddHandler("AliLoader", detName,
1015                                   loaderName, detName + "base",
1016                                   loaderName + "(const char*, TFolder*)");
1017         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1018       }
1019       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1020         fLoader[iDet] =
1021           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1022                                                  fRunLoader->GetEventFolder());
1023       }
1024       if (!fLoader[iDet]) {   // use default loader
1025         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1026       }
1027       if (!fLoader[iDet]) {
1028         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1029         if (fStopOnError) return NULL;
1030       } else {
1031         fRunLoader->AddLoader(fLoader[iDet]);
1032         fRunLoader->CdGAFile();
1033         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1034         fRunLoader->Write(0, TObject::kOverwrite);
1035       }
1036     }
1037   }
1038
1039   return reconstructor;
1040 }
1041
1042 //_____________________________________________________________________________
1043 Bool_t AliReconstruction::CreateVertexer()
1044 {
1045 // create the vertexer
1046
1047   fVertexer = NULL;
1048   AliReconstructor* itsReconstructor = GetReconstructor(0);
1049   if (itsReconstructor) {
1050     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1051   }
1052   if (!fVertexer) {
1053     AliWarning("couldn't create a vertexer for ITS");
1054     if (fStopOnError) return kFALSE;
1055   }
1056
1057   return kTRUE;
1058 }
1059
1060 //_____________________________________________________________________________
1061 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1062 {
1063 // create the trackers
1064
1065   TString detStr = detectors;
1066   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1067     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1068     AliReconstructor* reconstructor = GetReconstructor(iDet);
1069     if (!reconstructor) continue;
1070     TString detName = fgkDetectorName[iDet];
1071     if (detName == "HLT") {
1072       fRunHLTTracking = kTRUE;
1073       continue;
1074     }
1075
1076     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1077     if (!fTracker[iDet] && (iDet < 7)) {
1078       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1079       if (fStopOnError) return kFALSE;
1080     }
1081   }
1082
1083   return kTRUE;
1084 }
1085
1086 //_____________________________________________________________________________
1087 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1088 {
1089 // delete trackers and the run loader and close and delete the file
1090
1091   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1092     delete fReconstructor[iDet];
1093     fReconstructor[iDet] = NULL;
1094     fLoader[iDet] = NULL;
1095     delete fTracker[iDet];
1096     fTracker[iDet] = NULL;
1097   }
1098   delete fVertexer;
1099   fVertexer = NULL;
1100
1101   delete fRunLoader;
1102   fRunLoader = NULL;
1103   delete fRawReader;
1104   fRawReader = NULL;
1105
1106   if (file) {
1107     file->Close();
1108     delete file;
1109   }
1110
1111   if (fileOld) {
1112     fileOld->Close();
1113     delete fileOld;
1114     gSystem->Unlink("AliESDs.old.root");
1115   }
1116 }
1117
1118
1119 //_____________________________________________________________________________
1120 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1121 {
1122 // read the ESD event from a file
1123
1124   if (!esd) return kFALSE;
1125   char fileName[256];
1126   sprintf(fileName, "ESD_%d.%d_%s.root",
1127           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1128   if (gSystem->AccessPathName(fileName)) return kFALSE;
1129
1130   AliInfo(Form("reading ESD from file %s", fileName));
1131   AliDebug(1, Form("reading ESD from file %s", fileName));
1132   TFile* file = TFile::Open(fileName);
1133   if (!file || !file->IsOpen()) {
1134     AliError(Form("opening %s failed", fileName));
1135     delete file;
1136     return kFALSE;
1137   }
1138
1139   gROOT->cd();
1140   delete esd;
1141   esd = (AliESD*) file->Get("ESD");
1142   file->Close();
1143   delete file;
1144   return kTRUE;
1145 }
1146
1147 //_____________________________________________________________________________
1148 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1149 {
1150 // write the ESD event to a file
1151
1152   if (!esd) return;
1153   char fileName[256];
1154   sprintf(fileName, "ESD_%d.%d_%s.root",
1155           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1156
1157   AliDebug(1, Form("writing ESD to file %s", fileName));
1158   TFile* file = TFile::Open(fileName, "recreate");
1159   if (!file || !file->IsOpen()) {
1160     AliError(Form("opening %s failed", fileName));
1161   } else {
1162     esd->Write("ESD");
1163     file->Close();
1164   }
1165   delete file;
1166 }
1167
1168
1169
1170
1171 //_____________________________________________________________________________
1172 void AliReconstruction::CreateTag(TFile* file)
1173 {
1174   Int_t ntrack;
1175   Int_t NProtons, NKaons, NPions, NMuons, NElectrons;
1176   Int_t Npos, Nneg, Nneutr;
1177   Int_t NK0s, Nneutrons, Npi0s, Ngamas;
1178   Int_t Nch1GeV, Nch3GeV, Nch10GeV;
1179   Int_t Nmu1GeV, Nmu3GeV, Nmu10GeV;
1180   Int_t Nel1GeV, Nel3GeV, Nel10GeV;
1181   Float_t MaxPt = .0, MeanPt = .0, TotalP = .0;
1182
1183   AliRunTag *tag = new AliRunTag();
1184   AliDetectorTag *detTag = new AliDetectorTag();
1185   AliEventTag *evTag = new AliEventTag();
1186   TTree ttag("T","A Tree with event tags");
1187   TBranch * btag = ttag.Branch("AliTAG", "AliRunTag", &tag);
1188   btag->SetCompressionLevel(9);
1189
1190   AliInfo(Form("Creating the tags......."));
1191
1192   if (!file || !file->IsOpen()) {
1193     AliError(Form("opening failed"));
1194     delete file;
1195     return ;
1196   }
1197
1198   TTree *t = (TTree*) file->Get("esdTree");
1199   TBranch * b = t->GetBranch("ESD");
1200   AliESD *esd = 0;
1201   b->SetAddress(&esd);
1202
1203   tag->SetRunId(esd->GetRunNumber());
1204
1205   Int_t firstEvent = 0,lastEvent = 0;
1206   Int_t i_NumberOfEvents = b->GetEntries();
1207   for (Int_t i_EventNumber = 0; i_EventNumber < i_NumberOfEvents; i_EventNumber++)
1208     {
1209       ntrack = 0;
1210       Npos = 0;
1211       Nneg = 0;
1212       Nneutr =0;
1213       NK0s = 0;
1214       Nneutrons = 0;
1215       Npi0s = 0;
1216       Ngamas = 0;
1217       NProtons = 0;
1218       NKaons = 0;
1219       NPions = 0;
1220       NMuons = 0;
1221       NElectrons = 0;
1222       Nch1GeV = 0;
1223       Nch3GeV = 0;
1224       Nch10GeV = 0;
1225       Nmu1GeV = 0;
1226       Nmu3GeV = 0;
1227       Nmu10GeV = 0;
1228       Nel1GeV = 0;
1229       Nel3GeV = 0;
1230       Nel10GeV = 0;
1231       MaxPt = .0;
1232       MeanPt = .0;
1233       TotalP = .0;
1234
1235       b->GetEntry(i_EventNumber);
1236       const AliESDVertex * VertexIn = esd->GetVertex();
1237
1238       for (Int_t i_TrackNumber = 0; i_TrackNumber < esd->GetNumberOfTracks(); i_TrackNumber++)
1239         {
1240           AliESDtrack * ESDTrack = esd->GetTrack(i_TrackNumber);
1241           UInt_t status = ESDTrack->GetStatus();
1242
1243           //select only tracks with ITS refit
1244           if ((status&AliESDtrack::kITSrefit)==0) continue;
1245
1246           //select only tracks with TPC refit-->remove extremely high Pt tracks
1247           if ((status&AliESDtrack::kTPCrefit)==0) continue;
1248
1249           //select only tracks with the "combined PID"
1250           if ((status&AliESDtrack::kESDpid)==0) continue;
1251                   Double_t p[3];
1252           ESDTrack->GetPxPyPz(p);
1253           Double_t P = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1254           Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1255           TotalP += P;
1256           MeanPt += fPt;
1257           if(fPt > MaxPt)
1258             MaxPt = fPt;
1259
1260           if(ESDTrack->GetSign() > 0)
1261             {
1262               Npos++;
1263               if(fPt > 1.0)
1264                 Nch1GeV++;
1265               if(fPt > 3.0)
1266                 Nch3GeV++;
1267               if(fPt > 10.0)
1268                 Nch10GeV++;
1269             }
1270           if(ESDTrack->GetSign() < 0)
1271             {
1272               Nneg++;
1273               if(fPt > 1.0)
1274                 Nch1GeV++;
1275               if(fPt > 3.0)
1276                 Nch3GeV++;
1277               if(fPt > 10.0)
1278                 Nch10GeV++;
1279             }
1280           if(ESDTrack->GetSign() == 0)
1281             Nneutr++;
1282
1283           //PID
1284           Double_t prob[10];
1285           ESDTrack->GetESDpid(prob);
1286
1287           //K0s
1288           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]))
1289             NK0s++;
1290           //neutrons
1291           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]))
1292             Nneutrons++;
1293           //pi0s
1294           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]))
1295             Npi0s++;
1296           //gamas
1297           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]))
1298             Ngamas++;
1299           //protons
1300           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]))
1301             NProtons++;
1302           //kaons
1303           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]))
1304             NKaons++;
1305           //kaons
1306           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]))
1307             NPions++;
1308           //muons
1309           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]))
1310             {
1311               NMuons++;
1312               if(fPt > 1.0)
1313                 Nmu1GeV++;
1314               if(fPt > 3.0)
1315                 Nmu3GeV++;
1316               if(fPt > 10.0)
1317                 Nmu10GeV++;
1318             }
1319           //electrons
1320           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]))
1321             {
1322               NElectrons++;
1323               if(fPt > 1.0)
1324                 Nel1GeV++;
1325               if(fPt > 3.0)
1326                 Nel3GeV++;
1327               if(fPt > 10.0)
1328                 Nel10GeV++;
1329             }
1330
1331
1332
1333           ntrack++;
1334         }//track loop
1335       // Fill the event tags
1336       if(ntrack != 0)
1337         MeanPt = MeanPt/ntrack;
1338
1339       evTag->SetEventId(i_EventNumber+1);
1340       evTag->SetVertexX(VertexIn->GetXv());
1341       evTag->SetVertexY(VertexIn->GetYv());
1342       evTag->SetVertexZ(VertexIn->GetZv());
1343
1344       evTag->SetT0VertexZ(esd->GetT0zVertex());
1345
1346       evTag->SetTrigger(esd->GetTrigger());
1347
1348       evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
1349       evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
1350       evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1351       evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1352
1353
1354       evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1355       evTag->SetNumOfPosTracks(Npos);
1356       evTag->SetNumOfNegTracks(Nneg);
1357       evTag->SetNumOfNeutrTracks(Nneutr);
1358
1359       evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1360       evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1361       evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1362       evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1363
1364       evTag->SetNumOfProtons(NProtons);
1365       evTag->SetNumOfKaons(NKaons);
1366       evTag->SetNumOfPions(NPions);
1367       evTag->SetNumOfMuons(NMuons);
1368       evTag->SetNumOfElectrons(NElectrons);
1369       evTag->SetNumOfPhotons(Ngamas);
1370       evTag->SetNumOfPi0s(Npi0s);
1371       evTag->SetNumOfNeutrons(Nneutrons);
1372       evTag->SetNumOfKaon0s(NK0s);
1373
1374       evTag->SetNumOfChargedAbove1GeV(Nch1GeV);
1375       evTag->SetNumOfChargedAbove3GeV(Nch3GeV);
1376       evTag->SetNumOfChargedAbove10GeV(Nch10GeV);
1377       evTag->SetNumOfMuonsAbove1GeV(Nmu1GeV);
1378       evTag->SetNumOfMuonsAbove3GeV(Nmu3GeV);
1379       evTag->SetNumOfMuonsAbove10GeV(Nmu10GeV);
1380       evTag->SetNumOfElectronsAbove1GeV(Nel1GeV);
1381       evTag->SetNumOfElectronsAbove3GeV(Nel3GeV);
1382       evTag->SetNumOfElectronsAbove10GeV(Nel10GeV);
1383
1384       evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
1385       evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
1386
1387       evTag->SetTotalMomentum(TotalP);
1388       evTag->SetMeanPt(MeanPt);
1389       evTag->SetMaxPt(MaxPt);
1390
1391       tag->AddEventTag(*evTag);
1392     }
1393   lastEvent = i_NumberOfEvents;
1394
1395   ttag.Fill();
1396   tag->Clear();
1397
1398   char fileName[256];
1399   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1400           tag->GetRunId(),firstEvent,lastEvent );
1401   AliInfo(Form("writing tags to file %s", fileName));
1402   AliDebug(1, Form("writing tags to file %s", fileName));
1403
1404   TFile* ftag = TFile::Open(fileName, "recreate");
1405   ftag->cd();
1406   ttag.Write();
1407   ftag->Close();
1408   file->cd();
1409   delete tag;
1410   delete detTag;
1411   delete evTag;
1412 }
1413