STEER/AliReconstruction.cxx

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