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 // In case of raw-data reconstruction the user can modify the default        //
  51 // number of events per digits/clusters/tracks file. In case the option      //
  52 // is not used the number is set 1. In case the user provides 0, than        //
  53 // the number of events is equal to the number of events inside the          //
  54 // raw-data file (i.e. one digits/clusters/tracks file):                     //
  55 //                                                                           //
  56 //   rec.SetNumberOfEventsPerFile(...);                                      //
  57 //                                                                           //
  58 //                                                                           //
  59 // The name of the galice file can be changed from the default               //
  60 // "galice.root" by passing it as argument to the AliReconstruction          //
  61 // constructor or by                                                         //
  62 //                                                                           //
  63 //   rec.SetGAliceFile("...");                                               //
  64 //                                                                           //
  65 // The local reconstruction can be switched on or off for individual         //
  66 // detectors by                                                              //
  67 //                                                                           //
  68 //   rec.SetRunLocalReconstruction("...");                                   //
  69 //                                                                           //
  70 // The argument is a (case sensitive) string with the names of the           //
  71 // detectors separated by a space. The special string "ALL" selects all      //
  72 // available detectors. This is the default.                                 //
  73 //                                                                           //
  74 // The reconstruction of the primary vertex position can be switched off by  //
  75 //                                                                           //
  76 //   rec.SetRunVertexFinder(kFALSE);                                         //
  77 //                                                                           //
  78 // The tracking and the creation of ESD tracks can be switched on for        //
  79 // selected detectors by                                                     //
  80 //                                                                           //
  81 //   rec.SetRunTracking("...");                                              //
  82 //                                                                           //
  83 // Uniform/nonuniform field tracking switches (default: uniform field)       //
  84 //                                                                           //
  85 //   rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
  86 //                                                                           //
  87 // The filling of additional ESD information can be steered by               //
  88 //                                                                           //
  89 //   rec.SetFillESD("...");                                                  //
  90 //                                                                           //
  91 // Again, for both methods the string specifies the list of detectors.       //
  92 // The default is "ALL".                                                     //
  93 //                                                                           //
  94 // The call of the shortcut method                                           //
  95 //                                                                           //
  96 //   rec.SetRunReconstruction("...");                                        //
  97 //                                                                           //
  98 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
  99 // SetFillESD with the same detector selecting string as argument.           //
 100 //                                                                           //
 101 // The reconstruction requires digits or raw data as input. For the creation //
 102 // of digits and raw data have a look at the class AliSimulation.            //
 103 //                                                                           //
 104 // For debug purposes the method SetCheckPointLevel can be used. If the      //
 105 // argument is greater than 0, files with ESD events will be written after   //
 106 // selected steps of the reconstruction for each event:                      //
 107 //   level 1: after tracking and after filling of ESD (final)                //
 108 //   level 2: in addition after each tracking step                           //
 109 //   level 3: in addition after the filling of ESD for each detector         //
 110 // If a final check point file exists for an event, this event will be       //
 111 // skipped in the reconstruction. The tracking and the filling of ESD for    //
 112 // a detector will be skipped as well, if the corresponding check point      //
 113 // file exists. The ESD event will then be loaded from the file instead.     //
 114 //                                                                           //
 115 ///////////////////////////////////////////////////////////////////////////////
 116
 117 #include <TArrayF.h>
 118 #include <TFile.h>
 119 #include <TSystem.h>
 120 #include <TROOT.h>
 121 #include <TPluginManager.h>
 122 #include <TStopwatch.h>
 123 #include <TGeoManager.h>
 124 #include <TLorentzVector.h>
 125
 126 #include "AliReconstruction.h"
 127 #include "AliReconstructor.h"
 128 #include "AliLog.h"
 129 #include "AliRunLoader.h"
 130 #include "AliRun.h"
 131 #include "AliRawReaderFile.h"
 132 #include "AliRawReaderDate.h"
 133 #include "AliRawReaderRoot.h"
 134 #include "AliRawEventHeaderBase.h"
 135 #include "AliESD.h"
 136 #include "AliESDfriend.h"
 137 #include "AliESDVertex.h"
 138 #include "AliMultiplicity.h"
 139 #include "AliTracker.h"
 140 #include "AliVertexer.h"
 141 #include "AliVertexerTracks.h"
 142 #include "AliV0vertexer.h"
 143 #include "AliCascadeVertexer.h"
 144 #include "AliHeader.h"
 145 #include "AliGenEventHeader.h"
 146 #include "AliPID.h"
 147 #include "AliESDpid.h"
 148 #include "AliESDtrack.h"
 149
 150 #include "AliRunTag.h"
 151 #include "AliDetectorTag.h"
 152 #include "AliEventTag.h"
 153
 154 #include "AliTrackPointArray.h"
 155 #include "AliCDBManager.h"
 156 #include "AliCDBEntry.h"
 157 #include "AliAlignObj.h"
 158
 159 #include "AliCentralTrigger.h"
 160 #include "AliCTPRawStream.h"
 161
 162 ClassImp(AliReconstruction)
 163
 164
 165 //_____________________________________________________________________________
 166 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
 167
 168 //_____________________________________________________________________________
 169 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
 170                                      const char* name, const char* title) :
 171   TNamed(name, title),
 172
 173   fUniformField(kTRUE),
 174   fRunVertexFinder(kTRUE),
 175   fRunHLTTracking(kFALSE),
 176   fRunMuonTracking(kFALSE),
 177   fStopOnError(kFALSE),
 178   fWriteAlignmentData(kFALSE),
 179   fWriteESDfriend(kFALSE),
 180   fWriteAOD(kFALSE),
 181   fFillTriggerESD(kTRUE),
 182
 183   fRunLocalReconstruction("ALL"),
 184   fRunTracking("ALL"),
 185   fFillESD("ALL"),
 186   fGAliceFileName(gAliceFilename),
 187   fInput(""),
 188   fEquipIdMap(""),
 189   fFirstEvent(0),
 190   fLastEvent(-1),
 191   fNumberOfEventsPerFile(1),
 192   fCheckPointLevel(0),
 193   fOptions(),
 194   fLoadAlignFromCDB(kTRUE),
 195   fLoadAlignData("ALL"),
 196
 197   fRunLoader(NULL),
 198   fRawReader(NULL),
 199
 200   fVertexer(NULL),
 201   fDiamondProfile(NULL),
 202
 203   fAlignObjArray(NULL),
 204   fCDBUri(cdbUri),
 205   fSpecCDBUri()
 206 {
 207 // create reconstruction object with default parameters
 208
 209   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 210     fReconstructor[iDet] = NULL;
 211     fLoader[iDet] = NULL;
 212     fTracker[iDet] = NULL;
 213   }
 214   AliPID pid;
 215 }
 216
 217 //_____________________________________________________________________________
 218 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
 219   TNamed(rec),
 220
 221   fUniformField(rec.fUniformField),
 222   fRunVertexFinder(rec.fRunVertexFinder),
 223   fRunHLTTracking(rec.fRunHLTTracking),
 224   fRunMuonTracking(rec.fRunMuonTracking),
 225   fStopOnError(rec.fStopOnError),
 226   fWriteAlignmentData(rec.fWriteAlignmentData),
 227   fWriteESDfriend(rec.fWriteESDfriend),
 228   fWriteAOD(rec.fWriteAOD),
 229   fFillTriggerESD(rec.fFillTriggerESD),
 230
 231   fRunLocalReconstruction(rec.fRunLocalReconstruction),
 232   fRunTracking(rec.fRunTracking),
 233   fFillESD(rec.fFillESD),
 234   fGAliceFileName(rec.fGAliceFileName),
 235   fInput(rec.fInput),
 236   fEquipIdMap(rec.fEquipIdMap),
 237   fFirstEvent(rec.fFirstEvent),
 238   fLastEvent(rec.fLastEvent),
 239   fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
 240   fCheckPointLevel(0),
 241   fOptions(),
 242   fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
 243   fLoadAlignData(rec.fLoadAlignData),
 244
 245   fRunLoader(NULL),
 246   fRawReader(NULL),
 247
 248   fVertexer(NULL),
 249   fDiamondProfile(NULL),
 250
 251   fAlignObjArray(rec.fAlignObjArray),
 252   fCDBUri(rec.fCDBUri),
 253   fSpecCDBUri()
 254 {
 255 // copy constructor
 256
 257   for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
 258     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
 259   }
 260   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 261     fReconstructor[iDet] = NULL;
 262     fLoader[iDet] = NULL;
 263     fTracker[iDet] = NULL;
 264   }
 265   for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
 266     if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
 267   }
 268 }
 269
 270 //_____________________________________________________________________________
 271 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 272 {
 273 // assignment operator
 274
 275   this->~AliReconstruction();
 276   new(this) AliReconstruction(rec);
 277   return *this;
 278 }
 279
 280 //_____________________________________________________________________________
 281 AliReconstruction::~AliReconstruction()
 282 {
 283 // clean up
 284
 285   CleanUp();
 286   fOptions.Delete();
 287   fSpecCDBUri.Delete();
 288 }
 289
 290 //_____________________________________________________________________________
 291 void AliReconstruction::InitCDBStorage()
 292 {
 293 // activate a default CDB storage
 294 // First check if we have any CDB storage set, because it is used
 295 // to retrieve the calibration and alignment constants
 296
 297   AliCDBManager* man = AliCDBManager::Instance();
 298   if (man->IsDefaultStorageSet())
 299   {
 300     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 301     AliWarning("Default CDB storage has been already set !");
 302     AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
 303     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 304     fCDBUri = "";
 305   }
 306   else {
 307     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 308     AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
 309     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 310     man->SetDefaultStorage(fCDBUri);
 311   }
 312
 313   // Now activate the detector specific CDB storage locations
 314   for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
 315     TObject* obj = fSpecCDBUri[i];
 316     if (!obj) continue;
 317     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 318     AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
 319     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 320     man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
 321   }
 322   man->Print();
 323 }
 324
 325 //_____________________________________________________________________________
 326 void AliReconstruction::SetDefaultStorage(const char* uri) {
 327 // Store the desired default CDB storage location
 328 // Activate it later within the Run() method
 329
 330   fCDBUri = uri;
 331
 332 }
 333
 334 //_____________________________________________________________________________
 335 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
 336 // Store a detector-specific CDB storage location
 337 // Activate it later within the Run() method
 338
 339   AliCDBPath aPath(calibType);
 340   if(!aPath.IsValid()){
 341         // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
 342         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 343                 if(!strcmp(calibType, fgkDetectorName[iDet])) {
 344                         aPath.SetPath(Form("%s/*", calibType));
 345                         AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
 346                         break;
 347                 }
 348         }
 349         if(!aPath.IsValid()){
 350                 AliError(Form("Not a valid path or detector: %s", calibType));
 351                 return;
 352         }
 353   }
 354
 355   // check that calibType refers to a "valid" detector name
 356   Bool_t isDetector = kFALSE;
 357   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 358     TString detName = fgkDetectorName[iDet];
 359     if(aPath.GetLevel0() == detName) {
 360         isDetector = kTRUE;
 361         break;
 362     }
 363   }
 364
 365   if(!isDetector) {
 366         AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
 367         return;
 368   }
 369
 370   TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
 371   if (obj) fSpecCDBUri.Remove(obj);
 372   fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
 373
 374 }
 375
 376 //_____________________________________________________________________________
 377 Bool_t AliReconstruction::SetRunNumber()
 378 {
 379   // The method is called in Run() in order
 380   // to set a correct run number.
 381   // In case of raw data reconstruction the
 382   // run number is taken from the raw data header
 383
 384   if(AliCDBManager::Instance()->GetRun() < 0) {
 385     if (!fRunLoader) {
 386       AliError("No run loader is found !");
 387       return kFALSE;
 388     }
 389     // read run number from gAlice
 390     if(fRunLoader->GetAliRun())
 391       AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
 392     else {
 393       if(fRawReader) {
 394         if(fRawReader->NextEvent()) {
 395           AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
 396           fRawReader->RewindEvents();
 397         }
 398         else {
 399           AliError("No raw-data events found !");
 400           return kFALSE;
 401         }
 402       }
 403       else {
 404         AliError("Neither gAlice nor RawReader objects are found !");
 405         return kFALSE;
 406       }
 407     }
 408     AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
 409   }
 410   return kTRUE;
 411 }
 412
 413 //_____________________________________________________________________________
 414 Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
 415 {
 416   // Read collection of alignment objects (AliAlignObj derived) saved
 417   // in the TClonesArray ClArrayName and apply them to the geometry
 418   // manager singleton.
 419   //
 420   alObjArray->Sort();
 421   Int_t nvols = alObjArray->GetEntriesFast();
 422
 423   Bool_t flag = kTRUE;
 424
 425   for(Int_t j=0; j<nvols; j++)
 426     {
 427       AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
 428       if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
 429     }
 430
 431   if (AliDebugLevelClass() >= 1) {
 432     gGeoManager->GetTopNode()->CheckOverlaps(20);
 433     TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
 434     if(ovexlist->GetEntriesFast()){
 435       AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
 436    }
 437   }
 438
 439   return flag;
 440
 441 }
 442
 443 //_____________________________________________________________________________
 444 Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
 445 {
 446   // Fills array of single detector's alignable objects from CDB
 447
 448   AliDebug(2, Form("Loading alignment data for detector: %s",detName));
 449
 450   AliCDBEntry *entry;
 451
 452   AliCDBPath path(detName,"Align","Data");
 453
 454   entry=AliCDBManager::Instance()->Get(path.GetPath());
 455   if(!entry){
 456         AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
 457         return kFALSE;
 458   }
 459   entry->SetOwner(1);
 460   TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
 461   alignArray->SetOwner(0);
 462   AliDebug(2,Form("Found %d alignment objects for %s",
 463                         alignArray->GetEntries(),detName));
 464
 465   AliAlignObj *alignObj=0;
 466   TIter iter(alignArray);
 467
 468   // loop over align objects in detector
 469   while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
 470         fAlignObjArray->Add(alignObj);
 471   }
 472   // delete entry --- Don't delete, it is cached!
 473
 474   AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
 475   return kTRUE;
 476
 477 }
 478
 479 //_____________________________________________________________________________
 480 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
 481 {
 482   // Read the alignment objects from CDB.
 483   // Each detector is supposed to have the
 484   // alignment objects in DET/Align/Data CDB path.
 485   // All the detector objects are then collected,
 486   // sorted by geometry level (starting from ALIC) and
 487   // then applied to the TGeo geometry.
 488   // Finally an overlaps check is performed.
 489
 490   // Load alignment data from CDB and fill fAlignObjArray
 491   if(fLoadAlignFromCDB){
 492         if(!fAlignObjArray) fAlignObjArray = new TObjArray();
 493
 494         //fAlignObjArray->RemoveAll();
 495         fAlignObjArray->Clear();
 496         fAlignObjArray->SetOwner(0);
 497
 498         TString detStr = detectors;
 499         TString dataNotLoaded="";
 500         TString dataLoaded="";
 501
 502         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 503           if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 504           if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
 505             dataNotLoaded += fgkDetectorName[iDet];
 506             dataNotLoaded += " ";
 507           } else {
 508             dataLoaded += fgkDetectorName[iDet];
 509             dataLoaded += " ";
 510           }
 511         } // end loop over detectors
 512
 513         if ((detStr.CompareTo("ALL") == 0)) detStr = "";
 514         dataNotLoaded += detStr;
 515         if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
 516                           dataLoaded.Data()));
 517         if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
 518                           dataNotLoaded.Data()));
 519   } // fLoadAlignFromCDB flag
 520
 521   // Check if the array with alignment objects was
 522   // provided by the user. If yes, apply the objects
 523   // to the present TGeo geometry
 524   if (fAlignObjArray) {
 525     if (gGeoManager && gGeoManager->IsClosed()) {
 526       if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
 527         AliError("The misalignment of one or more volumes failed!"
 528                  "Compare the list of simulated detectors and the list of detector alignment data!");
 529         return kFALSE;
 530       }
 531     }
 532     else {
 533       AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
 534       return kFALSE;
 535     }
 536   }
 537
 538   delete fAlignObjArray; fAlignObjArray=0;
 539   return kTRUE;
 540 }
 541
 542 //_____________________________________________________________________________
 543 void AliReconstruction::SetGAliceFile(const char* fileName)
 544 {
 545 // set the name of the galice file
 546
 547   fGAliceFileName = fileName;
 548 }
 549
 550 //_____________________________________________________________________________
 551 void AliReconstruction::SetOption(const char* detector, const char* option)
 552 {
 553 // set options for the reconstruction of a detector
 554
 555   TObject* obj = fOptions.FindObject(detector);
 556   if (obj) fOptions.Remove(obj);
 557   fOptions.Add(new TNamed(detector, option));
 558 }
 559
 560
 561 //_____________________________________________________________________________
 562 Bool_t AliReconstruction::Run(const char* input)
 563 {
 564 // run the reconstruction
 565
 566   // set the input
 567   if (!input) input = fInput.Data();
 568   TString fileName(input);
 569   if (fileName.EndsWith("/")) {
 570     fRawReader = new AliRawReaderFile(fileName);
 571   } else if (fileName.EndsWith(".root")) {
 572     fRawReader = new AliRawReaderRoot(fileName);
 573   } else if (!fileName.IsNull()) {
 574     fRawReader = new AliRawReaderDate(fileName);
 575     fRawReader->SelectEvents(7);
 576   }
 577   if (!fEquipIdMap.IsNull() && fRawReader)
 578     fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
 579
 580
 581   // get the run loader
 582   if (!InitRunLoader()) return kFALSE;
 583
 584   // Initialize the CDB storage
 585   InitCDBStorage();
 586
 587   // Set run number in CDBManager (if it is not already set by the user)
 588   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
 589
 590   // Import ideal TGeo geometry and apply misalignment
 591   if (!gGeoManager) {
 592     TString geom(gSystem->DirName(fGAliceFileName));
 593     geom += "/geometry.root";
 594     TGeoManager::Import(geom.Data());
 595     if (!gGeoManager) if (fStopOnError) return kFALSE;
 596   }
 597
 598   AliCDBManager* man = AliCDBManager::Instance();
 599   if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
 600
 601   // local reconstruction
 602   if (!fRunLocalReconstruction.IsNull()) {
 603     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
 604       if (fStopOnError) {CleanUp(); return kFALSE;}
 605     }
 606   }
 607 //  if (!fRunVertexFinder && fRunTracking.IsNull() &&
 608 //      fFillESD.IsNull()) return kTRUE;
 609
 610   // get vertexer
 611   if (fRunVertexFinder && !CreateVertexer()) {
 612     if (fStopOnError) {
 613       CleanUp();
 614       return kFALSE;
 615     }
 616   }
 617
 618   // get trackers
 619   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
 620     if (fStopOnError) {
 621       CleanUp();
 622       return kFALSE;
 623     }
 624   }
 625
 626
 627   TStopwatch stopwatch;
 628   stopwatch.Start();
 629
 630   // get the possibly already existing ESD file and tree
 631   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
 632   TFile* fileOld = NULL;
 633   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
 634   if (!gSystem->AccessPathName("AliESDs.root")){
 635     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
 636     fileOld = TFile::Open("AliESDs.old.root");
 637     if (fileOld && fileOld->IsOpen()) {
 638       treeOld = (TTree*) fileOld->Get("esdTree");
 639       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
 640       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
 641       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
 642     }
 643   }
 644
 645   // create the ESD output file and tree
 646   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 647   if (!file->IsOpen()) {
 648     AliError("opening AliESDs.root failed");
 649     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 650   }
 651   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 652   tree->Branch("ESD", "AliESD", &esd);
 653   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
 654   hlttree->Branch("ESD", "AliESD", &hltesd);
 655   delete esd; delete hltesd;
 656   esd = NULL; hltesd = NULL;
 657
 658   // create the branch with ESD additions
 659   AliESDfriend *esdf=0;
 660   if (fWriteESDfriend) {
 661      TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
 662      br->SetFile("AliESDfriends.root");
 663   }
 664
 665   AliVertexerTracks tVertexer;
 666   if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
 667
 668   // loop over events
 669   if (fRawReader) fRawReader->RewindEvents();
 670
 671   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 672     if (fRawReader) fRawReader->NextEvent();
 673     if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
 674       // copy old ESD to the new one
 675       if (treeOld) {
 676         treeOld->SetBranchAddress("ESD", &esd);
 677         treeOld->GetEntry(iEvent);
 678       }
 679       tree->Fill();
 680       if (hlttreeOld) {
 681         hlttreeOld->SetBranchAddress("ESD", &hltesd);
 682         hlttreeOld->GetEntry(iEvent);
 683       }
 684       hlttree->Fill();
 685       continue;
 686     }
 687
 688     AliInfo(Form("processing event %d", iEvent));
 689     fRunLoader->GetEvent(iEvent);
 690
 691     char aFileName[256];
 692     sprintf(aFileName, "ESD_%d.%d_final.root",
 693             fRunLoader->GetHeader()->GetRun(),
 694             fRunLoader->GetHeader()->GetEventNrInRun());
 695     if (!gSystem->AccessPathName(aFileName)) continue;
 696
 697     // local reconstruction
 698     if (!fRunLocalReconstruction.IsNull()) {
 699       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
 700         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 701       }
 702     }
 703
 704     esd = new AliESD; hltesd = new AliESD;
 705     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 706     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 707     esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
 708     hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
 709
 710     // Set magnetic field from the tracker
 711     esd->SetMagneticField(AliTracker::GetBz());
 712     hltesd->SetMagneticField(AliTracker::GetBz());
 713
 714     // vertex finder
 715     if (fRunVertexFinder) {
 716       if (!ReadESD(esd, "vertex")) {
 717         if (!RunVertexFinder(esd)) {
 718           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 719         }
 720         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 721       }
 722     }
 723
 724     // HLT tracking
 725     if (!fRunTracking.IsNull()) {
 726       if (fRunHLTTracking) {
 727         hltesd->SetVertex(esd->GetVertex());
 728         if (!RunHLTTracking(hltesd)) {
 729           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 730         }
 731       }
 732     }
 733
 734     // Muon tracking
 735     if (!fRunTracking.IsNull()) {
 736       if (fRunMuonTracking) {
 737         if (!RunMuonTracking()) {
 738           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 739         }
 740       }
 741     }
 742
 743     // barrel tracking
 744     if (!fRunTracking.IsNull()) {
 745       if (!fRunMuonTracking) {
 746         if (!ReadESD(esd, "tracking")) {
 747           if (!RunTracking(esd)) {
 748             if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 749           }
 750           if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 751         }
 752       }
 753     }
 754     // fill ESD
 755     if (!fFillESD.IsNull()) {
 756       if (!FillESD(esd, fFillESD)) {
 757         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 758       }
 759     }
 760     // fill Event header information from the RawEventHeader
 761     if (fRawReader){FillRawEventHeaderESD(esd);}
 762
 763     // combined PID
 764     AliESDpid::MakePID(esd);
 765     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 766
 767     if (fFillTriggerESD) {
 768       if (!ReadESD(esd, "trigger")) {
 769         if (!FillTriggerESD(esd)) {
 770           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 771         }
 772         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
 773       }
 774     }
 775
 776     //Try to improve the reconstructed primary vertex position using the tracks
 777     AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
 778     if (pvtx)
 779     if (pvtx->GetStatus()) {
 780        // Store the improved primary vertex
 781        esd->SetPrimaryVertex(pvtx);
 782        // Propagate the tracks to the DCA to the improved primary vertex
 783        Double_t somethingbig = 777.;
 784        Double_t bz = esd->GetMagneticField();
 785        Int_t nt=esd->GetNumberOfTracks();
 786        while (nt--) {
 787          AliESDtrack *t = esd->GetTrack(nt);
 788          t->RelateToVertex(pvtx, bz, somethingbig);
 789        }
 790     }
 791
 792     {
 793     // V0 finding
 794     AliV0vertexer vtxer;
 795     vtxer.Tracks2V0vertices(esd);
 796
 797     // Cascade finding
 798     AliCascadeVertexer cvtxer;
 799     cvtxer.V0sTracks2CascadeVertices(esd);
 800     }
 801
 802     // write ESD
 803     if (fWriteESDfriend) {
 804        esdf=new AliESDfriend();
 805        esd->GetESDfriend(esdf);
 806     }
 807     tree->Fill();
 808
 809     // write HLT ESD
 810     hlttree->Fill();
 811
 812     if (fCheckPointLevel > 0)  WriteESD(esd, "final");
 813
 814     delete esd; delete esdf; delete hltesd;
 815     esd = NULL; esdf=NULL; hltesd = NULL;
 816   }
 817
 818   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
 819                stopwatch.RealTime(),stopwatch.CpuTime()));
 820
 821   file->cd();
 822   if (fWriteESDfriend)
 823     tree->SetBranchStatus("ESDfriend*",0);
 824   tree->Write();
 825   hlttree->Write();
 826
 827   if (fWriteAOD) {
 828     CreateAOD(file);
 829   }
 830
 831   // Create tags for the events in the ESD tree (the ESD tree is always present)
 832   // In case of empty events the tags will contain dummy values
 833   CreateTag(file);
 834   CleanUp(file, fileOld);
 835
 836   return kTRUE;
 837 }
 838
 839
 840 //_____________________________________________________________________________
 841 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 842 {
 843 // run the local reconstruction
 844
 845   TStopwatch stopwatch;
 846   stopwatch.Start();
 847
 848   AliCDBManager* man = AliCDBManager::Instance();
 849   Bool_t origCache = man->GetCacheFlag();
 850
 851   TString detStr = detectors;
 852   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 853     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 854     AliReconstructor* reconstructor = GetReconstructor(iDet);
 855     if (!reconstructor) continue;
 856     if (reconstructor->HasLocalReconstruction()) continue;
 857
 858     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 859     TStopwatch stopwatchDet;
 860     stopwatchDet.Start();
 861
 862     AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
 863
 864     man->SetCacheFlag(kTRUE);
 865     TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
 866     man->GetAll(calibPath); // entries are cached!
 867
 868     if (fRawReader) {
 869       fRawReader->RewindEvents();
 870       reconstructor->Reconstruct(fRunLoader, fRawReader);
 871     } else {
 872       reconstructor->Reconstruct(fRunLoader);
 873     }
 874     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 875                  fgkDetectorName[iDet],
 876                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 877
 878     // unload calibration data
 879     man->ClearCache();
 880   }
 881
 882   man->SetCacheFlag(origCache);
 883
 884   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 885     AliError(Form("the following detectors were not found: %s",
 886                   detStr.Data()));
 887     if (fStopOnError) return kFALSE;
 888   }
 889
 890   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 891                stopwatch.RealTime(),stopwatch.CpuTime()));
 892
 893   return kTRUE;
 894 }
 895
 896 //_____________________________________________________________________________
 897 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 898 {
 899 // run the local reconstruction
 900
 901   TStopwatch stopwatch;
 902   stopwatch.Start();
 903
 904   TString detStr = detectors;
 905   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 906     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 907     AliReconstructor* reconstructor = GetReconstructor(iDet);
 908     if (!reconstructor) continue;
 909     AliLoader* loader = fLoader[iDet];
 910
 911     // conversion of digits
 912     if (fRawReader && reconstructor->HasDigitConversion()) {
 913       AliInfo(Form("converting raw data digits into root objects for %s",
 914                    fgkDetectorName[iDet]));
 915       TStopwatch stopwatchDet;
 916       stopwatchDet.Start();
 917       loader->LoadDigits("update");
 918       loader->CleanDigits();
 919       loader->MakeDigitsContainer();
 920       TTree* digitsTree = loader->TreeD();
 921       reconstructor->ConvertDigits(fRawReader, digitsTree);
 922       loader->WriteDigits("OVERWRITE");
 923       loader->UnloadDigits();
 924       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 925                    fgkDetectorName[iDet],
 926                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 927     }
 928
 929     // local reconstruction
 930     if (!reconstructor->HasLocalReconstruction()) continue;
 931     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 932     TStopwatch stopwatchDet;
 933     stopwatchDet.Start();
 934     loader->LoadRecPoints("update");
 935     loader->CleanRecPoints();
 936     loader->MakeRecPointsContainer();
 937     TTree* clustersTree = loader->TreeR();
 938     if (fRawReader && !reconstructor->HasDigitConversion()) {
 939       reconstructor->Reconstruct(fRawReader, clustersTree);
 940     } else {
 941       loader->LoadDigits("read");
 942       TTree* digitsTree = loader->TreeD();
 943       if (!digitsTree) {
 944         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 945         if (fStopOnError) return kFALSE;
 946       } else {
 947         reconstructor->Reconstruct(digitsTree, clustersTree);
 948       }
 949       loader->UnloadDigits();
 950     }
 951     loader->WriteRecPoints("OVERWRITE");
 952     loader->UnloadRecPoints();
 953     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
 954                     fgkDetectorName[iDet],
 955                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 956   }
 957
 958   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 959     AliError(Form("the following detectors were not found: %s",
 960                   detStr.Data()));
 961     if (fStopOnError) return kFALSE;
 962   }
 963
 964   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 965                stopwatch.RealTime(),stopwatch.CpuTime()));
 966
 967   return kTRUE;
 968 }
 969
 970 //_____________________________________________________________________________
 971 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 972 {
 973 // run the barrel tracking
 974
 975   TStopwatch stopwatch;
 976   stopwatch.Start();
 977
 978   AliESDVertex* vertex = NULL;
 979   Double_t vtxPos[3] = {0, 0, 0};
 980   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 981   TArrayF mcVertex(3);
 982   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 983     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 984     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 985   }
 986
 987   if (fVertexer) {
 988     AliInfo("running the ITS vertex finder");
 989     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 990     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 991     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 992     if(!vertex){
 993       AliWarning("Vertex not found");
 994       vertex = new AliESDVertex();
 995       vertex->SetName("default");
 996     }
 997     else {
 998       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 999       vertex->SetName("reconstructed");
1000     }
1001
1002   } else {
1003     AliInfo("getting the primary vertex from MC");
1004     vertex = new AliESDVertex(vtxPos, vtxErr);
1005   }
1006
1007   if (vertex) {
1008     vertex->GetXYZ(vtxPos);
1009     vertex->GetSigmaXYZ(vtxErr);
1010   } else {
1011     AliWarning("no vertex reconstructed");
1012     vertex = new AliESDVertex(vtxPos, vtxErr);
1013   }
1014   esd->SetVertex(vertex);
1015   // if SPD multiplicity has been determined, it is stored in the ESD
1016   AliMultiplicity *mult= fVertexer->GetMultiplicity();
1017   if(mult)esd->SetMultiplicity(mult);
1018
1019   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1020     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1021   }
1022   delete vertex;
1023
1024   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1025                stopwatch.RealTime(),stopwatch.CpuTime()));
1026
1027   return kTRUE;
1028 }
1029
1030 //_____________________________________________________________________________
1031 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1032 {
1033 // run the HLT barrel tracking
1034
1035   TStopwatch stopwatch;
1036   stopwatch.Start();
1037
1038   if (!fRunLoader) {
1039     AliError("Missing runLoader!");
1040     return kFALSE;
1041   }
1042
1043   AliInfo("running HLT tracking");
1044
1045   // Get a pointer to the HLT reconstructor
1046   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1047   if (!reconstructor) return kFALSE;
1048
1049   // TPC + ITS
1050   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1051     TString detName = fgkDetectorName[iDet];
1052     AliDebug(1, Form("%s HLT tracking", detName.Data()));
1053     reconstructor->SetOption(detName.Data());
1054     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1055     if (!tracker) {
1056       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1057       if (fStopOnError) return kFALSE;
1058       continue;
1059     }
1060     Double_t vtxPos[3];
1061     Double_t vtxErr[3]={0.005,0.005,0.010};
1062     const AliESDVertex *vertex = esd->GetVertex();
1063     vertex->GetXYZ(vtxPos);
1064     tracker->SetVertex(vtxPos,vtxErr);
1065     if(iDet != 1) {
1066       fLoader[iDet]->LoadRecPoints("read");
1067       TTree* tree = fLoader[iDet]->TreeR();
1068       if (!tree) {
1069         AliError(Form("Can't get the %s cluster tree", detName.Data()));
1070         return kFALSE;
1071       }
1072       tracker->LoadClusters(tree);
1073     }
1074     if (tracker->Clusters2Tracks(esd) != 0) {
1075       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1076       return kFALSE;
1077     }
1078     if(iDet != 1) {
1079       tracker->UnloadClusters();
1080     }
1081     delete tracker;
1082   }
1083
1084   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1085                stopwatch.RealTime(),stopwatch.CpuTime()));
1086
1087   return kTRUE;
1088 }
1089
1090 //_____________________________________________________________________________
1091 Bool_t AliReconstruction::RunMuonTracking()
1092 {
1093 // run the muon spectrometer tracking
1094
1095   TStopwatch stopwatch;
1096   stopwatch.Start();
1097
1098   if (!fRunLoader) {
1099     AliError("Missing runLoader!");
1100     return kFALSE;
1101   }
1102   Int_t iDet = 7; // for MUON
1103
1104   AliInfo("is running...");
1105
1106   // Get a pointer to the MUON reconstructor
1107   AliReconstructor *reconstructor = GetReconstructor(iDet);
1108   if (!reconstructor) return kFALSE;
1109
1110
1111   TString detName = fgkDetectorName[iDet];
1112   AliDebug(1, Form("%s tracking", detName.Data()));
1113   AliTracker *tracker =  reconstructor->CreateTracker(fRunLoader);
1114   if (!tracker) {
1115     AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1116     return kFALSE;
1117   }
1118
1119   // create Tracks
1120   fLoader[iDet]->LoadTracks("update");
1121   fLoader[iDet]->CleanTracks();
1122   fLoader[iDet]->MakeTracksContainer();
1123
1124   // read RecPoints
1125   fLoader[iDet]->LoadRecPoints("read");
1126
1127   if (!tracker->Clusters2Tracks(0x0)) {
1128     AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1129     return kFALSE;
1130   }
1131   fLoader[iDet]->UnloadRecPoints();
1132
1133   fLoader[iDet]->WriteTracks("OVERWRITE");
1134   fLoader[iDet]->UnloadTracks();
1135
1136   delete tracker;
1137
1138
1139   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1140                stopwatch.RealTime(),stopwatch.CpuTime()));
1141
1142   return kTRUE;
1143 }
1144
1145
1146 //_____________________________________________________________________________
1147 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1148 {
1149 // run the barrel tracking
1150
1151   TStopwatch stopwatch;
1152   stopwatch.Start();
1153
1154   AliInfo("running tracking");
1155
1156   // pass 1: TPC + ITS inwards
1157   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1158     if (!fTracker[iDet]) continue;
1159     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1160
1161     // load clusters
1162     fLoader[iDet]->LoadRecPoints("read");
1163     TTree* tree = fLoader[iDet]->TreeR();
1164     if (!tree) {
1165       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1166       return kFALSE;
1167     }
1168     fTracker[iDet]->LoadClusters(tree);
1169
1170     // run tracking
1171     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1172       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1173       return kFALSE;
1174     }
1175     if (fCheckPointLevel > 1) {
1176       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1177     }
1178     // preliminary PID in TPC needed by the ITS tracker
1179     if (iDet == 1) {
1180       GetReconstructor(1)->FillESD(fRunLoader, esd);
1181       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1182       AliESDpid::MakePID(esd);
1183     }
1184   }
1185
1186   // pass 2: ALL backwards
1187   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1188     if (!fTracker[iDet]) continue;
1189     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1190
1191     // load clusters
1192     if (iDet > 1) {     // all except ITS, TPC
1193       TTree* tree = NULL;
1194       fLoader[iDet]->LoadRecPoints("read");
1195       tree = fLoader[iDet]->TreeR();
1196       if (!tree) {
1197         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1198         return kFALSE;
1199       }
1200       fTracker[iDet]->LoadClusters(tree);
1201     }
1202
1203     // run tracking
1204     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1205       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1206       return kFALSE;
1207     }
1208     if (fCheckPointLevel > 1) {
1209       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1210     }
1211
1212     // unload clusters
1213     if (iDet > 2) {     // all except ITS, TPC, TRD
1214       fTracker[iDet]->UnloadClusters();
1215       fLoader[iDet]->UnloadRecPoints();
1216     }
1217     // updated PID in TPC needed by the ITS tracker -MI
1218     if (iDet == 1) {
1219       GetReconstructor(1)->FillESD(fRunLoader, esd);
1220       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1221       AliESDpid::MakePID(esd);
1222     }
1223   }
1224
1225   // write space-points to the ESD in case alignment data output
1226   // is switched on
1227   if (fWriteAlignmentData)
1228     WriteAlignmentData(esd);
1229
1230   // pass 3: TRD + TPC + ITS refit inwards
1231   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1232     if (!fTracker[iDet]) continue;
1233     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1234
1235     // run tracking
1236     if (fTracker[iDet]->RefitInward(esd) != 0) {
1237       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1238       return kFALSE;
1239     }
1240     if (fCheckPointLevel > 1) {
1241       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1242     }
1243
1244     // unload clusters
1245     fTracker[iDet]->UnloadClusters();
1246     fLoader[iDet]->UnloadRecPoints();
1247   }
1248   //
1249   // Propagate track to the vertex - if not done by ITS
1250   //
1251   Int_t ntracks = esd->GetNumberOfTracks();
1252   for (Int_t itrack=0; itrack<ntracks; itrack++){
1253     const Double_t kRadius  = 3;   // beam pipe radius
1254     const Double_t kMaxStep = 5;   // max step
1255     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1256     Double_t       fieldZ   = AliTracker::GetBz();  //
1257     AliESDtrack * track = esd->GetTrack(itrack);
1258     if (!track) continue;
1259     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1260     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1261     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1262   }
1263
1264   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1265                stopwatch.RealTime(),stopwatch.CpuTime()));
1266
1267   return kTRUE;
1268 }
1269
1270 //_____________________________________________________________________________
1271 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1272 {
1273 // fill the event summary data
1274
1275   TStopwatch stopwatch;
1276   stopwatch.Start();
1277   AliInfo("filling ESD");
1278
1279   TString detStr = detectors;
1280   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1281     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1282     AliReconstructor* reconstructor = GetReconstructor(iDet);
1283     if (!reconstructor) continue;
1284
1285     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1286       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1287       TTree* clustersTree = NULL;
1288       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1289         fLoader[iDet]->LoadRecPoints("read");
1290         clustersTree = fLoader[iDet]->TreeR();
1291         if (!clustersTree) {
1292           AliError(Form("Can't get the %s clusters tree",
1293                         fgkDetectorName[iDet]));
1294           if (fStopOnError) return kFALSE;
1295         }
1296       }
1297       if (fRawReader && !reconstructor->HasDigitConversion()) {
1298         reconstructor->FillESD(fRawReader, clustersTree, esd);
1299       } else {
1300         TTree* digitsTree = NULL;
1301         if (fLoader[iDet]) {
1302           fLoader[iDet]->LoadDigits("read");
1303           digitsTree = fLoader[iDet]->TreeD();
1304           if (!digitsTree) {
1305             AliError(Form("Can't get the %s digits tree",
1306                           fgkDetectorName[iDet]));
1307             if (fStopOnError) return kFALSE;
1308           }
1309         }
1310         reconstructor->FillESD(digitsTree, clustersTree, esd);
1311         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1312       }
1313       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1314         fLoader[iDet]->UnloadRecPoints();
1315       }
1316
1317       if (fRawReader) {
1318         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1319       } else {
1320         reconstructor->FillESD(fRunLoader, esd);
1321       }
1322       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1323     }
1324   }
1325
1326   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1327     AliError(Form("the following detectors were not found: %s",
1328                   detStr.Data()));
1329     if (fStopOnError) return kFALSE;
1330   }
1331
1332   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1333                stopwatch.RealTime(),stopwatch.CpuTime()));
1334
1335   return kTRUE;
1336 }
1337
1338 //_____________________________________________________________________________
1339 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1340 {
1341   // Reads the trigger decision which is
1342   // stored in Trigger.root file and fills
1343   // the corresponding esd entries
1344
1345   AliInfo("Filling trigger information into the ESD");
1346
1347   if (fRawReader) {
1348     AliCTPRawStream input(fRawReader);
1349     if (!input.Next()) {
1350       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1351       return kFALSE;
1352     }
1353     esd->SetTriggerMask(input.GetClassMask());
1354     esd->SetTriggerCluster(input.GetClusterMask());
1355   }
1356   else {
1357     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1358     if (runloader) {
1359       if (!runloader->LoadTrigger()) {
1360         AliCentralTrigger *aCTP = runloader->GetTrigger();
1361         esd->SetTriggerMask(aCTP->GetClassMask());
1362         esd->SetTriggerCluster(aCTP->GetClusterMask());
1363       }
1364       else {
1365         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1366         return kFALSE;
1367       }
1368     }
1369     else {
1370       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1371       return kFALSE;
1372     }
1373   }
1374
1375   return kTRUE;
1376 }
1377
1378
1379
1380
1381
1382 //_____________________________________________________________________________
1383 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1384 {
1385   //
1386   // Filling information from RawReader Header
1387   //
1388
1389   AliInfo("Filling information from RawReader Header");
1390   esd->SetBunchCrossNumber(0);
1391   esd->SetOrbitNumber(0);
1392   esd->SetTimeStamp(0);
1393   esd->SetEventType(0);
1394   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1395   if (eventHeader){
1396     esd->SetBunchCrossNumber((eventHeader->GetP("Id")[0]));
1397     esd->SetOrbitNumber((eventHeader->GetP("Id")[1]));
1398     esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1399     esd->SetEventType((eventHeader->Get("Type")));
1400   }
1401
1402   return kTRUE;
1403 }
1404
1405
1406 //_____________________________________________________________________________
1407 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1408 {
1409 // check whether detName is contained in detectors
1410 // if yes, it is removed from detectors
1411
1412   // check if all detectors are selected
1413   if ((detectors.CompareTo("ALL") == 0) ||
1414       detectors.BeginsWith("ALL ") ||
1415       detectors.EndsWith(" ALL") ||
1416       detectors.Contains(" ALL ")) {
1417     detectors = "ALL";
1418     return kTRUE;
1419   }
1420
1421   // search for the given detector
1422   Bool_t result = kFALSE;
1423   if ((detectors.CompareTo(detName) == 0) ||
1424       detectors.BeginsWith(detName+" ") ||
1425       detectors.EndsWith(" "+detName) ||
1426       detectors.Contains(" "+detName+" ")) {
1427     detectors.ReplaceAll(detName, "");
1428     result = kTRUE;
1429   }
1430
1431   // clean up the detectors string
1432   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1433   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1434   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1435
1436   return result;
1437 }
1438
1439 //_____________________________________________________________________________
1440 Bool_t AliReconstruction::InitRunLoader()
1441 {
1442 // get or create the run loader
1443
1444   if (gAlice) delete gAlice;
1445   gAlice = NULL;
1446
1447   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1448     // load all base libraries to get the loader classes
1449     TString libs = gSystem->GetLibraries();
1450     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1451       TString detName = fgkDetectorName[iDet];
1452       if (detName == "HLT") continue;
1453       if (libs.Contains("lib" + detName + "base.so")) continue;
1454       gSystem->Load("lib" + detName + "base.so");
1455     }
1456     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1457     if (!fRunLoader) {
1458       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1459       CleanUp();
1460       return kFALSE;
1461     }
1462     fRunLoader->CdGAFile();
1463     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1464       if (fRunLoader->LoadgAlice() == 0) {
1465         gAlice = fRunLoader->GetAliRun();
1466         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1467       }
1468     }
1469     if (!gAlice && !fRawReader) {
1470       AliError(Form("no gAlice object found in file %s",
1471                     fGAliceFileName.Data()));
1472       CleanUp();
1473       return kFALSE;
1474     }
1475
1476   } else {               // galice.root does not exist
1477     if (!fRawReader) {
1478       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1479       CleanUp();
1480       return kFALSE;
1481     }
1482     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1483                                     AliConfig::GetDefaultEventFolderName(),
1484                                     "recreate");
1485     if (!fRunLoader) {
1486       AliError(Form("could not create run loader in file %s",
1487                     fGAliceFileName.Data()));
1488       CleanUp();
1489       return kFALSE;
1490     }
1491     fRunLoader->MakeTree("E");
1492     Int_t iEvent = 0;
1493     while (fRawReader->NextEvent()) {
1494       fRunLoader->SetEventNumber(iEvent);
1495       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1496                                      iEvent, iEvent);
1497       fRunLoader->MakeTree("H");
1498       fRunLoader->TreeE()->Fill();
1499       iEvent++;
1500     }
1501     fRawReader->RewindEvents();
1502     if (fNumberOfEventsPerFile > 0)
1503       fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1504     else
1505       fRunLoader->SetNumberOfEventsPerFile(iEvent);
1506     fRunLoader->WriteHeader("OVERWRITE");
1507     fRunLoader->CdGAFile();
1508     fRunLoader->Write(0, TObject::kOverwrite);
1509 //    AliTracker::SetFieldMap(???);
1510   }
1511
1512   return kTRUE;
1513 }
1514
1515 //_____________________________________________________________________________
1516 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1517 {
1518 // get the reconstructor object and the loader for a detector
1519
1520   if (fReconstructor[iDet]) return fReconstructor[iDet];
1521
1522   // load the reconstructor object
1523   TPluginManager* pluginManager = gROOT->GetPluginManager();
1524   TString detName = fgkDetectorName[iDet];
1525   TString recName = "Ali" + detName + "Reconstructor";
1526   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1527
1528   if (detName == "HLT") {
1529     if (!gROOT->GetClass("AliLevel3")) {
1530       gSystem->Load("libAliHLTSrc.so");
1531       gSystem->Load("libAliHLTMisc.so");
1532       gSystem->Load("libAliHLTHough.so");
1533       gSystem->Load("libAliHLTComp.so");
1534     }
1535   }
1536
1537   AliReconstructor* reconstructor = NULL;
1538   // first check if a plugin is defined for the reconstructor
1539   TPluginHandler* pluginHandler =
1540     pluginManager->FindHandler("AliReconstructor", detName);
1541   // if not, add a plugin for it
1542   if (!pluginHandler) {
1543     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1544     TString libs = gSystem->GetLibraries();
1545     if (libs.Contains("lib" + detName + "base.so") ||
1546         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1547       pluginManager->AddHandler("AliReconstructor", detName,
1548                                 recName, detName + "rec", recName + "()");
1549     } else {
1550       pluginManager->AddHandler("AliReconstructor", detName,
1551                                 recName, detName, recName + "()");
1552     }
1553     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1554   }
1555   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1556     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1557   }
1558   if (reconstructor) {
1559     TObject* obj = fOptions.FindObject(detName.Data());
1560     if (obj) reconstructor->SetOption(obj->GetTitle());
1561     reconstructor->Init(fRunLoader);
1562     fReconstructor[iDet] = reconstructor;
1563   }
1564
1565   // get or create the loader
1566   if (detName != "HLT") {
1567     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1568     if (!fLoader[iDet]) {
1569       AliConfig::Instance()
1570         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1571                                 detName, detName);
1572       // first check if a plugin is defined for the loader
1573       pluginHandler =
1574         pluginManager->FindHandler("AliLoader", detName);
1575       // if not, add a plugin for it
1576       if (!pluginHandler) {
1577         TString loaderName = "Ali" + detName + "Loader";
1578         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1579         pluginManager->AddHandler("AliLoader", detName,
1580                                   loaderName, detName + "base",
1581                                   loaderName + "(const char*, TFolder*)");
1582         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1583       }
1584       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1585         fLoader[iDet] =
1586           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1587                                                  fRunLoader->GetEventFolder());
1588       }
1589       if (!fLoader[iDet]) {   // use default loader
1590         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1591       }
1592       if (!fLoader[iDet]) {
1593         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1594         if (fStopOnError) return NULL;
1595       } else {
1596         fRunLoader->AddLoader(fLoader[iDet]);
1597         fRunLoader->CdGAFile();
1598         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1599         fRunLoader->Write(0, TObject::kOverwrite);
1600       }
1601     }
1602   }
1603
1604   return reconstructor;
1605 }
1606
1607 //_____________________________________________________________________________
1608 Bool_t AliReconstruction::CreateVertexer()
1609 {
1610 // create the vertexer
1611
1612   fVertexer = NULL;
1613   AliReconstructor* itsReconstructor = GetReconstructor(0);
1614   if (itsReconstructor) {
1615     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1616   }
1617   if (!fVertexer) {
1618     AliWarning("couldn't create a vertexer for ITS");
1619     if (fStopOnError) return kFALSE;
1620   }
1621
1622   return kTRUE;
1623 }
1624
1625 //_____________________________________________________________________________
1626 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1627 {
1628 // create the trackers
1629
1630   TString detStr = detectors;
1631   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1632     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1633     AliReconstructor* reconstructor = GetReconstructor(iDet);
1634     if (!reconstructor) continue;
1635     TString detName = fgkDetectorName[iDet];
1636     if (detName == "HLT") {
1637       fRunHLTTracking = kTRUE;
1638       continue;
1639     }
1640     if (detName == "MUON") {
1641       fRunMuonTracking = kTRUE;
1642       continue;
1643     }
1644
1645
1646     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1647     if (!fTracker[iDet] && (iDet < 7)) {
1648       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1649       if (fStopOnError) return kFALSE;
1650     }
1651   }
1652
1653   return kTRUE;
1654 }
1655
1656 //_____________________________________________________________________________
1657 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1658 {
1659 // delete trackers and the run loader and close and delete the file
1660
1661   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1662     delete fReconstructor[iDet];
1663     fReconstructor[iDet] = NULL;
1664     fLoader[iDet] = NULL;
1665     delete fTracker[iDet];
1666     fTracker[iDet] = NULL;
1667   }
1668   delete fVertexer;
1669   fVertexer = NULL;
1670   delete fDiamondProfile;
1671   fDiamondProfile = NULL;
1672
1673   delete fRunLoader;
1674   fRunLoader = NULL;
1675   delete fRawReader;
1676   fRawReader = NULL;
1677
1678   if (file) {
1679     file->Close();
1680     delete file;
1681   }
1682
1683   if (fileOld) {
1684     fileOld->Close();
1685     delete fileOld;
1686     gSystem->Unlink("AliESDs.old.root");
1687   }
1688 }
1689
1690
1691 //_____________________________________________________________________________
1692 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1693 {
1694 // read the ESD event from a file
1695
1696   if (!esd) return kFALSE;
1697   char fileName[256];
1698   sprintf(fileName, "ESD_%d.%d_%s.root",
1699           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1700   if (gSystem->AccessPathName(fileName)) return kFALSE;
1701
1702   AliInfo(Form("reading ESD from file %s", fileName));
1703   AliDebug(1, Form("reading ESD from file %s", fileName));
1704   TFile* file = TFile::Open(fileName);
1705   if (!file || !file->IsOpen()) {
1706     AliError(Form("opening %s failed", fileName));
1707     delete file;
1708     return kFALSE;
1709   }
1710
1711   gROOT->cd();
1712   delete esd;
1713   esd = (AliESD*) file->Get("ESD");
1714   file->Close();
1715   delete file;
1716   return kTRUE;
1717 }
1718
1719 //_____________________________________________________________________________
1720 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1721 {
1722 // write the ESD event to a file
1723
1724   if (!esd) return;
1725   char fileName[256];
1726   sprintf(fileName, "ESD_%d.%d_%s.root",
1727           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1728
1729   AliDebug(1, Form("writing ESD to file %s", fileName));
1730   TFile* file = TFile::Open(fileName, "recreate");
1731   if (!file || !file->IsOpen()) {
1732     AliError(Form("opening %s failed", fileName));
1733   } else {
1734     esd->Write("ESD");
1735     file->Close();
1736   }
1737   delete file;
1738 }
1739
1740
1741
1742
1743 //_____________________________________________________________________________
1744 void AliReconstruction::CreateTag(TFile* file)
1745 {
1746   /////////////
1747   //muon code//
1748   ////////////
1749   Double_t fMUONMASS = 0.105658369;
1750   //Variables
1751   Double_t fX,fY,fZ ;
1752   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1753   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1754   Int_t fCharge;
1755   TLorentzVector fEPvector;
1756
1757   Float_t fZVertexCut = 10.0;
1758   Float_t fRhoVertexCut = 2.0;
1759
1760   Float_t fLowPtCut = 1.0;
1761   Float_t fHighPtCut = 3.0;
1762   Float_t fVeryHighPtCut = 10.0;
1763   ////////////
1764
1765   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1766
1767   // Creates the tags for all the events in a given ESD file
1768   Int_t ntrack;
1769   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1770   Int_t nPos, nNeg, nNeutr;
1771   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1772   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1773   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1774   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1775   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1776   Int_t fVertexflag;
1777   Int_t iRunNumber = 0;
1778   TString fVertexName("default");
1779
1780   AliRunTag *tag = new AliRunTag();
1781   AliEventTag *evTag = new AliEventTag();
1782   TTree ttag("T","A Tree with event tags");
1783   TBranch * btag = ttag.Branch("AliTAG", &tag);
1784   btag->SetCompressionLevel(9);
1785
1786   AliInfo(Form("Creating the tags......."));
1787
1788   if (!file || !file->IsOpen()) {
1789     AliError(Form("opening failed"));
1790     delete file;
1791     return ;
1792   }
1793   Int_t lastEvent = 0;
1794   TTree *t = (TTree*) file->Get("esdTree");
1795   TBranch * b = t->GetBranch("ESD");
1796   AliESD *esd = 0;
1797   b->SetAddress(&esd);
1798
1799   b->GetEntry(fFirstEvent);
1800   Int_t iInitRunNumber = esd->GetRunNumber();
1801
1802   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1803   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1804   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1805     ntrack = 0;
1806     nPos = 0;
1807     nNeg = 0;
1808     nNeutr =0;
1809     nK0s = 0;
1810     nNeutrons = 0;
1811     nPi0s = 0;
1812     nGamas = 0;
1813     nProtons = 0;
1814     nKaons = 0;
1815     nPions = 0;
1816     nMuons = 0;
1817     nElectrons = 0;
1818     nCh1GeV = 0;
1819     nCh3GeV = 0;
1820     nCh10GeV = 0;
1821     nMu1GeV = 0;
1822     nMu3GeV = 0;
1823     nMu10GeV = 0;
1824     nEl1GeV = 0;
1825     nEl3GeV = 0;
1826     nEl10GeV = 0;
1827     maxPt = .0;
1828     meanPt = .0;
1829     totalP = .0;
1830     fVertexflag = 0;
1831
1832     b->GetEntry(iEventNumber);
1833     iRunNumber = esd->GetRunNumber();
1834     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1835     const AliESDVertex * vertexIn = esd->GetVertex();
1836     if (!vertexIn) AliError("ESD has not defined vertex.");
1837     if (vertexIn) fVertexName = vertexIn->GetName();
1838     if(fVertexName != "default") fVertexflag = 1;
1839     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1840       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1841       UInt_t status = esdTrack->GetStatus();
1842
1843       //select only tracks with ITS refit
1844       if ((status&AliESDtrack::kITSrefit)==0) continue;
1845       //select only tracks with TPC refit
1846       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1847
1848       //select only tracks with the "combined PID"
1849       if ((status&AliESDtrack::kESDpid)==0) continue;
1850       Double_t p[3];
1851       esdTrack->GetPxPyPz(p);
1852       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1853       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1854       totalP += momentum;
1855       meanPt += fPt;
1856       if(fPt > maxPt) maxPt = fPt;
1857
1858       if(esdTrack->GetSign() > 0) {
1859         nPos++;
1860         if(fPt > fLowPtCut) nCh1GeV++;
1861         if(fPt > fHighPtCut) nCh3GeV++;
1862         if(fPt > fVeryHighPtCut) nCh10GeV++;
1863       }
1864       if(esdTrack->GetSign() < 0) {
1865         nNeg++;
1866         if(fPt > fLowPtCut) nCh1GeV++;
1867         if(fPt > fHighPtCut) nCh3GeV++;
1868         if(fPt > fVeryHighPtCut) nCh10GeV++;
1869       }
1870       if(esdTrack->GetSign() == 0) nNeutr++;
1871
1872       //PID
1873       Double_t prob[5];
1874       esdTrack->GetESDpid(prob);
1875
1876       Double_t rcc = 0.0;
1877       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1878       if(rcc == 0.0) continue;
1879       //Bayes' formula
1880       Double_t w[5];
1881       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1882
1883       //protons
1884       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1885       //kaons
1886       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1887       //pions
1888       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1889       //electrons
1890       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1891         nElectrons++;
1892         if(fPt > fLowPtCut) nEl1GeV++;
1893         if(fPt > fHighPtCut) nEl3GeV++;
1894         if(fPt > fVeryHighPtCut) nEl10GeV++;
1895       }
1896       ntrack++;
1897     }//track loop
1898
1899     /////////////
1900     //muon code//
1901     ////////////
1902     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1903     // loop over all reconstructed tracks (also first track of combination)
1904     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1905       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1906       if (muonTrack == 0x0) continue;
1907
1908       // Coordinates at vertex
1909       fZ = muonTrack->GetZ();
1910       fY = muonTrack->GetBendingCoor();
1911       fX = muonTrack->GetNonBendingCoor();
1912
1913       fThetaX = muonTrack->GetThetaX();
1914       fThetaY = muonTrack->GetThetaY();
1915
1916       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1917       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1918       fPxRec = fPzRec * TMath::Tan(fThetaX);
1919       fPyRec = fPzRec * TMath::Tan(fThetaY);
1920       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1921
1922       //ChiSquare of the track if needed
1923       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1924       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1925       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1926
1927       // total number of muons inside a vertex cut
1928       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1929         nMuons++;
1930         if(fEPvector.Pt() > fLowPtCut) {
1931           nMu1GeV++;
1932           if(fEPvector.Pt() > fHighPtCut) {
1933             nMu3GeV++;
1934             if (fEPvector.Pt() > fVeryHighPtCut) {
1935               nMu10GeV++;
1936             }
1937           }
1938         }
1939       }
1940     }//muon track loop
1941
1942     // Fill the event tags
1943     if(ntrack != 0)
1944       meanPt = meanPt/ntrack;
1945
1946     evTag->SetEventId(iEventNumber+1);
1947     if (vertexIn) {
1948       evTag->SetVertexX(vertexIn->GetXv());
1949       evTag->SetVertexY(vertexIn->GetYv());
1950       evTag->SetVertexZ(vertexIn->GetZv());
1951       evTag->SetVertexZError(vertexIn->GetZRes());
1952     }
1953     evTag->SetVertexFlag(fVertexflag);
1954
1955     evTag->SetT0VertexZ(esd->GetT0zVertex());
1956
1957     evTag->SetTriggerMask(esd->GetTriggerMask());
1958     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1959
1960     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1961     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1962     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1963     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1964     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1965     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1966
1967
1968     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1969     evTag->SetNumOfPosTracks(nPos);
1970     evTag->SetNumOfNegTracks(nNeg);
1971     evTag->SetNumOfNeutrTracks(nNeutr);
1972
1973     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1974     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1975     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1976     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1977
1978     evTag->SetNumOfProtons(nProtons);
1979     evTag->SetNumOfKaons(nKaons);
1980     evTag->SetNumOfPions(nPions);
1981     evTag->SetNumOfMuons(nMuons);
1982     evTag->SetNumOfElectrons(nElectrons);
1983     evTag->SetNumOfPhotons(nGamas);
1984     evTag->SetNumOfPi0s(nPi0s);
1985     evTag->SetNumOfNeutrons(nNeutrons);
1986     evTag->SetNumOfKaon0s(nK0s);
1987
1988     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1989     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1990     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1991     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1992     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1993     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1994     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1995     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1996     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1997
1998     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1999     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2000
2001     evTag->SetTotalMomentum(totalP);
2002     evTag->SetMeanPt(meanPt);
2003     evTag->SetMaxPt(maxPt);
2004
2005     tag->SetRunId(iInitRunNumber);
2006     tag->AddEventTag(*evTag);
2007   }
2008   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2009   else lastEvent = fLastEvent;
2010
2011   ttag.Fill();
2012   tag->Clear();
2013
2014   char fileName[256];
2015   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2016           tag->GetRunId(),fFirstEvent,lastEvent );
2017   AliInfo(Form("writing tags to file %s", fileName));
2018   AliDebug(1, Form("writing tags to file %s", fileName));
2019
2020   TFile* ftag = TFile::Open(fileName, "recreate");
2021   ftag->cd();
2022   ttag.Write();
2023   ftag->Close();
2024   file->cd();
2025   delete tag;
2026   delete evTag;
2027 }
2028
2029 //_____________________________________________________________________________
2030 void AliReconstruction::CreateAOD(TFile* esdFile)
2031 {
2032   // do nothing for now
2033
2034   return;
2035 }
2036
2037
2038 void AliReconstruction::WriteAlignmentData(AliESD* esd)
2039 {
2040   // Write space-points which are then used in the alignment procedures
2041   // For the moment only ITS, TRD and TPC
2042
2043   // Load TOF clusters
2044   if (fTracker[3]){
2045     fLoader[3]->LoadRecPoints("read");
2046     TTree* tree = fLoader[3]->TreeR();
2047     if (!tree) {
2048       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2049       return;
2050     }
2051     fTracker[3]->LoadClusters(tree);
2052   }
2053   Int_t ntracks = esd->GetNumberOfTracks();
2054   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2055     {
2056       AliESDtrack *track = esd->GetTrack(itrack);
2057       Int_t nsp = 0;
2058       Int_t idx[200];
2059       for (Int_t iDet = 3; iDet >= 0; iDet--)
2060         nsp += track->GetNcls(iDet);
2061       if (nsp) {
2062         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2063         track->SetTrackPointArray(sp);
2064         Int_t isptrack = 0;
2065         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2066           AliTracker *tracker = fTracker[iDet];
2067           if (!tracker) continue;
2068           Int_t nspdet = track->GetNcls(iDet);
2069           if (nspdet <= 0) continue;
2070           track->GetClusters(iDet,idx);
2071           AliTrackPoint p;
2072           Int_t isp = 0;
2073           Int_t isp2 = 0;
2074           while (isp < nspdet) {
2075             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2076             const Int_t kNTPCmax = 159;
2077             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2078             if (!isvalid) continue;
2079             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2080           }
2081         }
2082       }
2083     }
2084   if (fTracker[3]){
2085     fTracker[3]->UnloadClusters();
2086     fLoader[3]->UnloadRecPoints();
2087   }
2088 }