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 "AliGeomManager.h"
 155 #include "AliTrackPointArray.h"
 156 #include "AliCDBManager.h"
 157 #include "AliCDBEntry.h"
 158 #include "AliAlignObj.h"
 159
 160 #include "AliCentralTrigger.h"
 161 #include "AliCTPRawStream.h"
 162
 163 #include "AliAODEvent.h"
 164 #include "AliAODHeader.h"
 165 #include "AliAODTrack.h"
 166 #include "AliAODVertex.h"
 167 #include "AliAODCluster.h"
 168
 169
 170 ClassImp(AliReconstruction)
 171
 172
 173 //_____________________________________________________________________________
 174 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
 175
 176 //_____________________________________________________________________________
 177 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
 178                                      const char* name, const char* title) :
 179   TNamed(name, title),
 180
 181   fUniformField(kTRUE),
 182   fRunVertexFinder(kTRUE),
 183   fRunHLTTracking(kFALSE),
 184   fRunMuonTracking(kFALSE),
 185   fStopOnError(kFALSE),
 186   fWriteAlignmentData(kFALSE),
 187   fWriteESDfriend(kFALSE),
 188   fWriteAOD(kFALSE),
 189   fFillTriggerESD(kTRUE),
 190
 191   fRunLocalReconstruction("ALL"),
 192   fRunTracking("ALL"),
 193   fFillESD("ALL"),
 194   fGAliceFileName(gAliceFilename),
 195   fInput(""),
 196   fEquipIdMap(""),
 197   fFirstEvent(0),
 198   fLastEvent(-1),
 199   fNumberOfEventsPerFile(1),
 200   fCheckPointLevel(0),
 201   fOptions(),
 202   fLoadAlignFromCDB(kTRUE),
 203   fLoadAlignData("ALL"),
 204
 205   fRunLoader(NULL),
 206   fRawReader(NULL),
 207
 208   fVertexer(NULL),
 209   fDiamondProfile(NULL),
 210
 211   fAlignObjArray(NULL),
 212   fCDBUri(cdbUri),
 213   fSpecCDBUri()
 214 {
 215 // create reconstruction object with default parameters
 216
 217   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 218     fReconstructor[iDet] = NULL;
 219     fLoader[iDet] = NULL;
 220     fTracker[iDet] = NULL;
 221   }
 222   AliPID pid;
 223 }
 224
 225 //_____________________________________________________________________________
 226 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
 227   TNamed(rec),
 228
 229   fUniformField(rec.fUniformField),
 230   fRunVertexFinder(rec.fRunVertexFinder),
 231   fRunHLTTracking(rec.fRunHLTTracking),
 232   fRunMuonTracking(rec.fRunMuonTracking),
 233   fStopOnError(rec.fStopOnError),
 234   fWriteAlignmentData(rec.fWriteAlignmentData),
 235   fWriteESDfriend(rec.fWriteESDfriend),
 236   fWriteAOD(rec.fWriteAOD),
 237   fFillTriggerESD(rec.fFillTriggerESD),
 238
 239   fRunLocalReconstruction(rec.fRunLocalReconstruction),
 240   fRunTracking(rec.fRunTracking),
 241   fFillESD(rec.fFillESD),
 242   fGAliceFileName(rec.fGAliceFileName),
 243   fInput(rec.fInput),
 244   fEquipIdMap(rec.fEquipIdMap),
 245   fFirstEvent(rec.fFirstEvent),
 246   fLastEvent(rec.fLastEvent),
 247   fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
 248   fCheckPointLevel(0),
 249   fOptions(),
 250   fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
 251   fLoadAlignData(rec.fLoadAlignData),
 252
 253   fRunLoader(NULL),
 254   fRawReader(NULL),
 255
 256   fVertexer(NULL),
 257   fDiamondProfile(NULL),
 258
 259   fAlignObjArray(rec.fAlignObjArray),
 260   fCDBUri(rec.fCDBUri),
 261   fSpecCDBUri()
 262 {
 263 // copy constructor
 264
 265   for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
 266     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
 267   }
 268   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 269     fReconstructor[iDet] = NULL;
 270     fLoader[iDet] = NULL;
 271     fTracker[iDet] = NULL;
 272   }
 273   for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
 274     if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
 275   }
 276 }
 277
 278 //_____________________________________________________________________________
 279 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 280 {
 281 // assignment operator
 282
 283   this->~AliReconstruction();
 284   new(this) AliReconstruction(rec);
 285   return *this;
 286 }
 287
 288 //_____________________________________________________________________________
 289 AliReconstruction::~AliReconstruction()
 290 {
 291 // clean up
 292
 293   CleanUp();
 294   fOptions.Delete();
 295   fSpecCDBUri.Delete();
 296 }
 297
 298 //_____________________________________________________________________________
 299 void AliReconstruction::InitCDBStorage()
 300 {
 301 // activate a default CDB storage
 302 // First check if we have any CDB storage set, because it is used
 303 // to retrieve the calibration and alignment constants
 304
 305   AliCDBManager* man = AliCDBManager::Instance();
 306   if (man->IsDefaultStorageSet())
 307   {
 308     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 309     AliWarning("Default CDB storage has been already set !");
 310     AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
 311     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 312     fCDBUri = "";
 313   }
 314   else {
 315     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 316     AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
 317     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 318     man->SetDefaultStorage(fCDBUri);
 319   }
 320
 321   // Now activate the detector specific CDB storage locations
 322   for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
 323     TObject* obj = fSpecCDBUri[i];
 324     if (!obj) continue;
 325     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 326     AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
 327     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 328     man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
 329   }
 330   man->Print();
 331 }
 332
 333 //_____________________________________________________________________________
 334 void AliReconstruction::SetDefaultStorage(const char* uri) {
 335 // Store the desired default CDB storage location
 336 // Activate it later within the Run() method
 337
 338   fCDBUri = uri;
 339
 340 }
 341
 342 //_____________________________________________________________________________
 343 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
 344 // Store a detector-specific CDB storage location
 345 // Activate it later within the Run() method
 346
 347   AliCDBPath aPath(calibType);
 348   if(!aPath.IsValid()){
 349         // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
 350         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 351                 if(!strcmp(calibType, fgkDetectorName[iDet])) {
 352                         aPath.SetPath(Form("%s/*", calibType));
 353                         AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
 354                         break;
 355                 }
 356         }
 357         if(!aPath.IsValid()){
 358                 AliError(Form("Not a valid path or detector: %s", calibType));
 359                 return;
 360         }
 361   }
 362
 363   // check that calibType refers to a "valid" detector name
 364   Bool_t isDetector = kFALSE;
 365   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 366     TString detName = fgkDetectorName[iDet];
 367     if(aPath.GetLevel0() == detName) {
 368         isDetector = kTRUE;
 369         break;
 370     }
 371   }
 372
 373   if(!isDetector) {
 374         AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
 375         return;
 376   }
 377
 378   TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
 379   if (obj) fSpecCDBUri.Remove(obj);
 380   fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
 381
 382 }
 383
 384 //_____________________________________________________________________________
 385 Bool_t AliReconstruction::SetRunNumber()
 386 {
 387   // The method is called in Run() in order
 388   // to set a correct run number.
 389   // In case of raw data reconstruction the
 390   // run number is taken from the raw data header
 391
 392   if(AliCDBManager::Instance()->GetRun() < 0) {
 393     if (!fRunLoader) {
 394       AliError("No run loader is found !");
 395       return kFALSE;
 396     }
 397     // read run number from gAlice
 398     if(fRunLoader->GetAliRun())
 399       AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
 400     else {
 401       if(fRawReader) {
 402         if(fRawReader->NextEvent()) {
 403           AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
 404           fRawReader->RewindEvents();
 405         }
 406         else {
 407           AliError("No raw-data events found !");
 408           return kFALSE;
 409         }
 410       }
 411       else {
 412         AliError("Neither gAlice nor RawReader objects are found !");
 413         return kFALSE;
 414       }
 415     }
 416     AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
 417   }
 418   return kTRUE;
 419 }
 420
 421 //_____________________________________________________________________________
 422 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
 423 {
 424   // Read the alignment objects from CDB.
 425   // Each detector is supposed to have the
 426   // alignment objects in DET/Align/Data CDB path.
 427   // All the detector objects are then collected,
 428   // sorted by geometry level (starting from ALIC) and
 429   // then applied to the TGeo geometry.
 430   // Finally an overlaps check is performed.
 431
 432   // Load alignment data from CDB and fill fAlignObjArray
 433   if(fLoadAlignFromCDB){
 434
 435     TString detStr = detectors;
 436     TString loadAlObjsListOfDets = "";
 437
 438     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 439       if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 440       loadAlObjsListOfDets += fgkDetectorName[iDet];
 441       loadAlObjsListOfDets += " ";
 442     } // end loop over detectors
 443     (AliGeomManager::Instance())->ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
 444   }else{
 445     // Check if the array with alignment objects was
 446     // provided by the user. If yes, apply the objects
 447     // to the present TGeo geometry
 448     if (fAlignObjArray) {
 449       if (gGeoManager && gGeoManager->IsClosed()) {
 450         if ((AliGeomManager::Instance())->ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
 451           AliError("The misalignment of one or more volumes failed!"
 452                    "Compare the list of simulated detectors and the list of detector alignment data!");
 453           return kFALSE;
 454         }
 455       }
 456       else {
 457         AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
 458         return kFALSE;
 459       }
 460     }
 461   }
 462
 463   delete fAlignObjArray; fAlignObjArray=0;
 464
 465   // Update the TGeoPhysicalNodes
 466   gGeoManager->RefreshPhysicalNodes();
 467
 468   return kTRUE;
 469 }
 470
 471 //_____________________________________________________________________________
 472 void AliReconstruction::SetGAliceFile(const char* fileName)
 473 {
 474 // set the name of the galice file
 475
 476   fGAliceFileName = fileName;
 477 }
 478
 479 //_____________________________________________________________________________
 480 void AliReconstruction::SetOption(const char* detector, const char* option)
 481 {
 482 // set options for the reconstruction of a detector
 483
 484   TObject* obj = fOptions.FindObject(detector);
 485   if (obj) fOptions.Remove(obj);
 486   fOptions.Add(new TNamed(detector, option));
 487 }
 488
 489
 490 //_____________________________________________________________________________
 491 Bool_t AliReconstruction::Run(const char* input)
 492 {
 493 // run the reconstruction
 494
 495   // set the input
 496   if (!input) input = fInput.Data();
 497   TString fileName(input);
 498   if (fileName.EndsWith("/")) {
 499     fRawReader = new AliRawReaderFile(fileName);
 500   } else if (fileName.EndsWith(".root")) {
 501     fRawReader = new AliRawReaderRoot(fileName);
 502   } else if (!fileName.IsNull()) {
 503     fRawReader = new AliRawReaderDate(fileName);
 504     fRawReader->SelectEvents(7);
 505   }
 506   if (!fEquipIdMap.IsNull() && fRawReader)
 507     fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
 508
 509
 510   // get the run loader
 511   if (!InitRunLoader()) return kFALSE;
 512
 513   // Initialize the CDB storage
 514   InitCDBStorage();
 515
 516   // Set run number in CDBManager (if it is not already set by the user)
 517   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
 518
 519   // Import ideal TGeo geometry and apply misalignment
 520   if (!gGeoManager) {
 521     TString geom(gSystem->DirName(fGAliceFileName));
 522     geom += "/geometry.root";
 523     TGeoManager::Import(geom.Data());
 524     if (!gGeoManager) if (fStopOnError) return kFALSE;
 525   }
 526
 527   AliCDBManager* man = AliCDBManager::Instance();
 528   if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
 529
 530   // local reconstruction
 531   if (!fRunLocalReconstruction.IsNull()) {
 532     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
 533       if (fStopOnError) {CleanUp(); return kFALSE;}
 534     }
 535   }
 536 //  if (!fRunVertexFinder && fRunTracking.IsNull() &&
 537 //      fFillESD.IsNull()) return kTRUE;
 538
 539   // get vertexer
 540   if (fRunVertexFinder && !CreateVertexer()) {
 541     if (fStopOnError) {
 542       CleanUp();
 543       return kFALSE;
 544     }
 545   }
 546
 547   // get trackers
 548   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
 549     if (fStopOnError) {
 550       CleanUp();
 551       return kFALSE;
 552     }
 553   }
 554
 555
 556   TStopwatch stopwatch;
 557   stopwatch.Start();
 558
 559   // get the possibly already existing ESD file and tree
 560   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
 561   TFile* fileOld = NULL;
 562   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
 563   if (!gSystem->AccessPathName("AliESDs.root")){
 564     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
 565     fileOld = TFile::Open("AliESDs.old.root");
 566     if (fileOld && fileOld->IsOpen()) {
 567       treeOld = (TTree*) fileOld->Get("esdTree");
 568       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
 569       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
 570       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
 571     }
 572   }
 573
 574   // create the ESD output file and tree
 575   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 576   if (!file->IsOpen()) {
 577     AliError("opening AliESDs.root failed");
 578     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 579   }
 580   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 581   tree->Branch("ESD", "AliESD", &esd);
 582   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
 583   hlttree->Branch("ESD", "AliESD", &hltesd);
 584   delete esd; delete hltesd;
 585   esd = NULL; hltesd = NULL;
 586
 587   // create the branch with ESD additions
 588   AliESDfriend *esdf=0;
 589   if (fWriteESDfriend) {
 590      TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
 591      br->SetFile("AliESDfriends.root");
 592   }
 593
 594   // Get the diamond profile from OCDB
 595   AliCDBEntry* entry = AliCDBManager::Instance()
 596         ->Get("GRP/Calib/MeanVertex");
 597
 598   if(entry) {
 599         fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
 600   } else {
 601         AliError("No diamond profile found in OCDB!");
 602   }
 603
 604   AliVertexerTracks tVertexer(AliTracker::GetBz());
 605   if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
 606
 607   // loop over events
 608   if (fRawReader) fRawReader->RewindEvents();
 609
 610   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 611     if (fRawReader) fRawReader->NextEvent();
 612     if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
 613       // copy old ESD to the new one
 614       if (treeOld) {
 615         treeOld->SetBranchAddress("ESD", &esd);
 616         treeOld->GetEntry(iEvent);
 617       }
 618       tree->Fill();
 619       if (hlttreeOld) {
 620         hlttreeOld->SetBranchAddress("ESD", &hltesd);
 621         hlttreeOld->GetEntry(iEvent);
 622       }
 623       hlttree->Fill();
 624       continue;
 625     }
 626
 627     AliInfo(Form("processing event %d", iEvent));
 628     fRunLoader->GetEvent(iEvent);
 629
 630     char aFileName[256];
 631     sprintf(aFileName, "ESD_%d.%d_final.root",
 632             fRunLoader->GetHeader()->GetRun(),
 633             fRunLoader->GetHeader()->GetEventNrInRun());
 634     if (!gSystem->AccessPathName(aFileName)) continue;
 635
 636     // local reconstruction
 637     if (!fRunLocalReconstruction.IsNull()) {
 638       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
 639         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 640       }
 641     }
 642
 643     esd = new AliESD; hltesd = new AliESD;
 644     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 645     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 646     esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
 647     hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
 648
 649     // Set magnetic field from the tracker
 650     esd->SetMagneticField(AliTracker::GetBz());
 651     hltesd->SetMagneticField(AliTracker::GetBz());
 652
 653     // Fill raw-data error log into the ESD
 654     if (fRawReader) FillRawDataErrorLog(iEvent,esd);
 655
 656     // vertex finder
 657     if (fRunVertexFinder) {
 658       if (!ReadESD(esd, "vertex")) {
 659         if (!RunVertexFinder(esd)) {
 660           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 661         }
 662         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 663       }
 664     }
 665
 666     // HLT tracking
 667     if (!fRunTracking.IsNull()) {
 668       if (fRunHLTTracking) {
 669         hltesd->SetVertex(esd->GetVertex());
 670         if (!RunHLTTracking(hltesd)) {
 671           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 672         }
 673       }
 674     }
 675
 676     // Muon tracking
 677     if (!fRunTracking.IsNull()) {
 678       if (fRunMuonTracking) {
 679         if (!RunMuonTracking()) {
 680           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 681         }
 682       }
 683     }
 684
 685     // barrel tracking
 686     if (!fRunTracking.IsNull()) {
 687       if (!ReadESD(esd, "tracking")) {
 688         if (!RunTracking(esd)) {
 689           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 690         }
 691         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 692       }
 693     }
 694
 695     // fill ESD
 696     if (!fFillESD.IsNull()) {
 697       if (!FillESD(esd, fFillESD)) {
 698         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 699       }
 700     }
 701     // fill Event header information from the RawEventHeader
 702     if (fRawReader){FillRawEventHeaderESD(esd);}
 703
 704     // combined PID
 705     AliESDpid::MakePID(esd);
 706     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 707
 708     if (fFillTriggerESD) {
 709       if (!ReadESD(esd, "trigger")) {
 710         if (!FillTriggerESD(esd)) {
 711           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 712         }
 713         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
 714       }
 715     }
 716
 717     //Try to improve the reconstructed primary vertex position using the tracks
 718     AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
 719     if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
 720
 721     if (pvtx)
 722     if (pvtx->GetStatus()) {
 723        // Store the improved primary vertex
 724        esd->SetPrimaryVertex(pvtx);
 725        // Propagate the tracks to the DCA to the improved primary vertex
 726        Double_t somethingbig = 777.;
 727        Double_t bz = esd->GetMagneticField();
 728        Int_t nt=esd->GetNumberOfTracks();
 729        while (nt--) {
 730          AliESDtrack *t = esd->GetTrack(nt);
 731          t->RelateToVertex(pvtx, bz, somethingbig);
 732        }
 733     }
 734
 735     {
 736     // V0 finding
 737     AliV0vertexer vtxer;
 738     vtxer.Tracks2V0vertices(esd);
 739
 740     // Cascade finding
 741     AliCascadeVertexer cvtxer;
 742     cvtxer.V0sTracks2CascadeVertices(esd);
 743     }
 744
 745     // write ESD
 746     if (fWriteESDfriend) {
 747        esdf=new AliESDfriend();
 748        esd->GetESDfriend(esdf);
 749     }
 750     tree->Fill();
 751
 752     // write HLT ESD
 753     hlttree->Fill();
 754
 755     if (fCheckPointLevel > 0)  WriteESD(esd, "final");
 756
 757     delete esd; delete esdf; delete hltesd;
 758     esd = NULL; esdf=NULL; hltesd = NULL;
 759   }
 760
 761   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
 762                stopwatch.RealTime(),stopwatch.CpuTime()));
 763
 764   file->cd();
 765   if (fWriteESDfriend)
 766     tree->SetBranchStatus("ESDfriend*",0);
 767   tree->Write();
 768   hlttree->Write();
 769
 770   if (fWriteAOD) {
 771     TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
 772     ESDFile2AODFile(file, aodFile);
 773     aodFile->Close();
 774   }
 775
 776   // Create tags for the events in the ESD tree (the ESD tree is always present)
 777   // In case of empty events the tags will contain dummy values
 778   CreateTag(file);
 779   CleanUp(file, fileOld);
 780
 781   return kTRUE;
 782 }
 783
 784
 785 //_____________________________________________________________________________
 786 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 787 {
 788 // run the local reconstruction
 789
 790   TStopwatch stopwatch;
 791   stopwatch.Start();
 792
 793   AliCDBManager* man = AliCDBManager::Instance();
 794   Bool_t origCache = man->GetCacheFlag();
 795
 796   TString detStr = detectors;
 797   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 798     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 799     AliReconstructor* reconstructor = GetReconstructor(iDet);
 800     if (!reconstructor) continue;
 801     if (reconstructor->HasLocalReconstruction()) continue;
 802
 803     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 804     TStopwatch stopwatchDet;
 805     stopwatchDet.Start();
 806
 807     AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
 808
 809     man->SetCacheFlag(kTRUE);
 810     TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
 811     man->GetAll(calibPath); // entries are cached!
 812
 813     if (fRawReader) {
 814       fRawReader->RewindEvents();
 815       reconstructor->Reconstruct(fRunLoader, fRawReader);
 816     } else {
 817       reconstructor->Reconstruct(fRunLoader);
 818     }
 819     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 820                  fgkDetectorName[iDet],
 821                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 822
 823     // unload calibration data
 824     man->ClearCache();
 825   }
 826
 827   man->SetCacheFlag(origCache);
 828
 829   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 830     AliError(Form("the following detectors were not found: %s",
 831                   detStr.Data()));
 832     if (fStopOnError) return kFALSE;
 833   }
 834
 835   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 836                stopwatch.RealTime(),stopwatch.CpuTime()));
 837
 838   return kTRUE;
 839 }
 840
 841 //_____________________________________________________________________________
 842 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 843 {
 844 // run the local reconstruction
 845
 846   TStopwatch stopwatch;
 847   stopwatch.Start();
 848
 849   TString detStr = detectors;
 850   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 851     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 852     AliReconstructor* reconstructor = GetReconstructor(iDet);
 853     if (!reconstructor) continue;
 854     AliLoader* loader = fLoader[iDet];
 855
 856     // conversion of digits
 857     if (fRawReader && reconstructor->HasDigitConversion()) {
 858       AliInfo(Form("converting raw data digits into root objects for %s",
 859                    fgkDetectorName[iDet]));
 860       TStopwatch stopwatchDet;
 861       stopwatchDet.Start();
 862       loader->LoadDigits("update");
 863       loader->CleanDigits();
 864       loader->MakeDigitsContainer();
 865       TTree* digitsTree = loader->TreeD();
 866       reconstructor->ConvertDigits(fRawReader, digitsTree);
 867       loader->WriteDigits("OVERWRITE");
 868       loader->UnloadDigits();
 869       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 870                    fgkDetectorName[iDet],
 871                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 872     }
 873
 874     // local reconstruction
 875     if (!reconstructor->HasLocalReconstruction()) continue;
 876     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 877     TStopwatch stopwatchDet;
 878     stopwatchDet.Start();
 879     loader->LoadRecPoints("update");
 880     loader->CleanRecPoints();
 881     loader->MakeRecPointsContainer();
 882     TTree* clustersTree = loader->TreeR();
 883     if (fRawReader && !reconstructor->HasDigitConversion()) {
 884       reconstructor->Reconstruct(fRawReader, clustersTree);
 885     } else {
 886       loader->LoadDigits("read");
 887       TTree* digitsTree = loader->TreeD();
 888       if (!digitsTree) {
 889         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 890         if (fStopOnError) return kFALSE;
 891       } else {
 892         reconstructor->Reconstruct(digitsTree, clustersTree);
 893       }
 894       loader->UnloadDigits();
 895     }
 896     loader->WriteRecPoints("OVERWRITE");
 897     loader->UnloadRecPoints();
 898     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 899                  fgkDetectorName[iDet],
 900                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 901   }
 902
 903   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 904     AliError(Form("the following detectors were not found: %s",
 905                   detStr.Data()));
 906     if (fStopOnError) return kFALSE;
 907   }
 908
 909   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 910                stopwatch.RealTime(),stopwatch.CpuTime()));
 911
 912   return kTRUE;
 913 }
 914
 915 //_____________________________________________________________________________
 916 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 917 {
 918 // run the barrel tracking
 919
 920   TStopwatch stopwatch;
 921   stopwatch.Start();
 922
 923   AliESDVertex* vertex = NULL;
 924   Double_t vtxPos[3] = {0, 0, 0};
 925   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 926   TArrayF mcVertex(3);
 927   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 928     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 929     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 930   }
 931
 932   if (fVertexer) {
 933     if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
 934     AliInfo("running the ITS vertex finder");
 935     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 936     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 937     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 938     if(!vertex){
 939       AliWarning("Vertex not found");
 940       vertex = new AliESDVertex();
 941       vertex->SetName("default");
 942     }
 943     else {
 944       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 945       vertex->SetName("reconstructed");
 946     }
 947
 948   } else {
 949     AliInfo("getting the primary vertex from MC");
 950     vertex = new AliESDVertex(vtxPos, vtxErr);
 951   }
 952
 953   if (vertex) {
 954     vertex->GetXYZ(vtxPos);
 955     vertex->GetSigmaXYZ(vtxErr);
 956   } else {
 957     AliWarning("no vertex reconstructed");
 958     vertex = new AliESDVertex(vtxPos, vtxErr);
 959   }
 960   esd->SetVertex(vertex);
 961   // if SPD multiplicity has been determined, it is stored in the ESD
 962   AliMultiplicity *mult = fVertexer->GetMultiplicity();
 963   if(mult)esd->SetMultiplicity(mult);
 964
 965   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 966     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
 967   }
 968   delete vertex;
 969
 970   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 971                stopwatch.RealTime(),stopwatch.CpuTime()));
 972
 973   return kTRUE;
 974 }
 975
 976 //_____________________________________________________________________________
 977 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
 978 {
 979 // run the HLT barrel tracking
 980
 981   TStopwatch stopwatch;
 982   stopwatch.Start();
 983
 984   if (!fRunLoader) {
 985     AliError("Missing runLoader!");
 986     return kFALSE;
 987   }
 988
 989   AliInfo("running HLT tracking");
 990
 991   // Get a pointer to the HLT reconstructor
 992   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
 993   if (!reconstructor) return kFALSE;
 994
 995   // TPC + ITS
 996   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 997     TString detName = fgkDetectorName[iDet];
 998     AliDebug(1, Form("%s HLT tracking", detName.Data()));
 999     reconstructor->SetOption(detName.Data());
1000     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1001     if (!tracker) {
1002       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1003       if (fStopOnError) return kFALSE;
1004       continue;
1005     }
1006     Double_t vtxPos[3];
1007     Double_t vtxErr[3]={0.005,0.005,0.010};
1008     const AliESDVertex *vertex = esd->GetVertex();
1009     vertex->GetXYZ(vtxPos);
1010     tracker->SetVertex(vtxPos,vtxErr);
1011     if(iDet != 1) {
1012       fLoader[iDet]->LoadRecPoints("read");
1013       TTree* tree = fLoader[iDet]->TreeR();
1014       if (!tree) {
1015         AliError(Form("Can't get the %s cluster tree", detName.Data()));
1016         return kFALSE;
1017       }
1018       tracker->LoadClusters(tree);
1019     }
1020     if (tracker->Clusters2Tracks(esd) != 0) {
1021       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1022       return kFALSE;
1023     }
1024     if(iDet != 1) {
1025       tracker->UnloadClusters();
1026     }
1027     delete tracker;
1028   }
1029
1030   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1031                stopwatch.RealTime(),stopwatch.CpuTime()));
1032
1033   return kTRUE;
1034 }
1035
1036 //_____________________________________________________________________________
1037 Bool_t AliReconstruction::RunMuonTracking()
1038 {
1039 // run the muon spectrometer tracking
1040
1041   TStopwatch stopwatch;
1042   stopwatch.Start();
1043
1044   if (!fRunLoader) {
1045     AliError("Missing runLoader!");
1046     return kFALSE;
1047   }
1048   Int_t iDet = 7; // for MUON
1049
1050   AliInfo("is running...");
1051
1052   // Get a pointer to the MUON reconstructor
1053   AliReconstructor *reconstructor = GetReconstructor(iDet);
1054   if (!reconstructor) return kFALSE;
1055
1056
1057   TString detName = fgkDetectorName[iDet];
1058   AliDebug(1, Form("%s tracking", detName.Data()));
1059   AliTracker *tracker =  reconstructor->CreateTracker(fRunLoader);
1060   if (!tracker) {
1061     AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1062     return kFALSE;
1063   }
1064
1065   // create Tracks
1066   fLoader[iDet]->LoadTracks("update");
1067   fLoader[iDet]->CleanTracks();
1068   fLoader[iDet]->MakeTracksContainer();
1069
1070   // read RecPoints
1071   fLoader[iDet]->LoadRecPoints("read");
1072
1073   if (!tracker->Clusters2Tracks(0x0)) {
1074     AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1075     return kFALSE;
1076   }
1077   fLoader[iDet]->UnloadRecPoints();
1078
1079   fLoader[iDet]->WriteTracks("OVERWRITE");
1080   fLoader[iDet]->UnloadTracks();
1081
1082   delete tracker;
1083
1084
1085   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1086                stopwatch.RealTime(),stopwatch.CpuTime()));
1087
1088   return kTRUE;
1089 }
1090
1091
1092 //_____________________________________________________________________________
1093 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1094 {
1095 // run the barrel tracking
1096
1097   TStopwatch stopwatch;
1098   stopwatch.Start();
1099
1100   AliInfo("running tracking");
1101
1102   //Fill the ESD with the T0 info (will be used by the TOF)
1103   if (fReconstructor[11])
1104       GetReconstructor(11)->FillESD(fRunLoader, esd);
1105
1106   // pass 1: TPC + ITS inwards
1107   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1108     if (!fTracker[iDet]) continue;
1109     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1110
1111     // load clusters
1112     fLoader[iDet]->LoadRecPoints("read");
1113     TTree* tree = fLoader[iDet]->TreeR();
1114     if (!tree) {
1115       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1116       return kFALSE;
1117     }
1118     fTracker[iDet]->LoadClusters(tree);
1119
1120     // run tracking
1121     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1122       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1123       return kFALSE;
1124     }
1125     if (fCheckPointLevel > 1) {
1126       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1127     }
1128     // preliminary PID in TPC needed by the ITS tracker
1129     if (iDet == 1) {
1130       GetReconstructor(1)->FillESD(fRunLoader, esd);
1131       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1132       AliESDpid::MakePID(esd);
1133     }
1134   }
1135
1136   // pass 2: ALL backwards
1137   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1138     if (!fTracker[iDet]) continue;
1139     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1140
1141     // load clusters
1142     if (iDet > 1) {     // all except ITS, TPC
1143       TTree* tree = NULL;
1144       fLoader[iDet]->LoadRecPoints("read");
1145       tree = fLoader[iDet]->TreeR();
1146       if (!tree) {
1147         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1148         return kFALSE;
1149       }
1150       fTracker[iDet]->LoadClusters(tree);
1151     }
1152
1153     // run tracking
1154     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1155       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1156       return kFALSE;
1157     }
1158     if (fCheckPointLevel > 1) {
1159       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1160     }
1161
1162     // unload clusters
1163     if (iDet > 2) {     // all except ITS, TPC, TRD
1164       fTracker[iDet]->UnloadClusters();
1165       fLoader[iDet]->UnloadRecPoints();
1166     }
1167     // updated PID in TPC needed by the ITS tracker -MI
1168     if (iDet == 1) {
1169       GetReconstructor(1)->FillESD(fRunLoader, esd);
1170       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1171       AliESDpid::MakePID(esd);
1172     }
1173   }
1174
1175   // write space-points to the ESD in case alignment data output
1176   // is switched on
1177   if (fWriteAlignmentData)
1178     WriteAlignmentData(esd);
1179
1180   // pass 3: TRD + TPC + ITS refit inwards
1181   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1182     if (!fTracker[iDet]) continue;
1183     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1184
1185     // run tracking
1186     if (fTracker[iDet]->RefitInward(esd) != 0) {
1187       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1188       return kFALSE;
1189     }
1190     if (fCheckPointLevel > 1) {
1191       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1192     }
1193
1194     // unload clusters
1195     fTracker[iDet]->UnloadClusters();
1196     fLoader[iDet]->UnloadRecPoints();
1197   }
1198   //
1199   // Propagate track to the vertex - if not done by ITS
1200   //
1201   Int_t ntracks = esd->GetNumberOfTracks();
1202   for (Int_t itrack=0; itrack<ntracks; itrack++){
1203     const Double_t kRadius  = 3;   // beam pipe radius
1204     const Double_t kMaxStep = 5;   // max step
1205     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1206     Double_t       fieldZ   = AliTracker::GetBz();  //
1207     AliESDtrack * track = esd->GetTrack(itrack);
1208     if (!track) continue;
1209     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1210     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1211     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1212   }
1213
1214   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1215                stopwatch.RealTime(),stopwatch.CpuTime()));
1216
1217   return kTRUE;
1218 }
1219
1220 //_____________________________________________________________________________
1221 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1222 {
1223 // fill the event summary data
1224
1225   TStopwatch stopwatch;
1226   stopwatch.Start();
1227   AliInfo("filling ESD");
1228
1229   TString detStr = detectors;
1230   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1231     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1232     AliReconstructor* reconstructor = GetReconstructor(iDet);
1233     if (!reconstructor) continue;
1234
1235     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1236       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1237       TTree* clustersTree = NULL;
1238       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1239         fLoader[iDet]->LoadRecPoints("read");
1240         clustersTree = fLoader[iDet]->TreeR();
1241         if (!clustersTree) {
1242           AliError(Form("Can't get the %s clusters tree",
1243                         fgkDetectorName[iDet]));
1244           if (fStopOnError) return kFALSE;
1245         }
1246       }
1247       if (fRawReader && !reconstructor->HasDigitConversion()) {
1248         reconstructor->FillESD(fRawReader, clustersTree, esd);
1249       } else {
1250         TTree* digitsTree = NULL;
1251         if (fLoader[iDet]) {
1252           fLoader[iDet]->LoadDigits("read");
1253           digitsTree = fLoader[iDet]->TreeD();
1254           if (!digitsTree) {
1255             AliError(Form("Can't get the %s digits tree",
1256                           fgkDetectorName[iDet]));
1257             if (fStopOnError) return kFALSE;
1258           }
1259         }
1260         reconstructor->FillESD(digitsTree, clustersTree, esd);
1261         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1262       }
1263       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1264         fLoader[iDet]->UnloadRecPoints();
1265       }
1266
1267       if (fRawReader) {
1268         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1269       } else {
1270         reconstructor->FillESD(fRunLoader, esd);
1271       }
1272       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1273     }
1274   }
1275
1276   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1277     AliError(Form("the following detectors were not found: %s",
1278                   detStr.Data()));
1279     if (fStopOnError) return kFALSE;
1280   }
1281
1282   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1283                stopwatch.RealTime(),stopwatch.CpuTime()));
1284
1285   return kTRUE;
1286 }
1287
1288 //_____________________________________________________________________________
1289 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1290 {
1291   // Reads the trigger decision which is
1292   // stored in Trigger.root file and fills
1293   // the corresponding esd entries
1294
1295   AliInfo("Filling trigger information into the ESD");
1296
1297   if (fRawReader) {
1298     AliCTPRawStream input(fRawReader);
1299     if (!input.Next()) {
1300       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1301       return kFALSE;
1302     }
1303     esd->SetTriggerMask(input.GetClassMask());
1304     esd->SetTriggerCluster(input.GetClusterMask());
1305   }
1306   else {
1307     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1308     if (runloader) {
1309       if (!runloader->LoadTrigger()) {
1310         AliCentralTrigger *aCTP = runloader->GetTrigger();
1311         esd->SetTriggerMask(aCTP->GetClassMask());
1312         esd->SetTriggerCluster(aCTP->GetClusterMask());
1313       }
1314       else {
1315         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1316         return kFALSE;
1317       }
1318     }
1319     else {
1320       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1321       return kFALSE;
1322     }
1323   }
1324
1325   return kTRUE;
1326 }
1327
1328
1329
1330
1331
1332 //_____________________________________________________________________________
1333 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1334 {
1335   //
1336   // Filling information from RawReader Header
1337   //
1338
1339   AliInfo("Filling information from RawReader Header");
1340   esd->SetBunchCrossNumber(0);
1341   esd->SetOrbitNumber(0);
1342   esd->SetPeriodNumber(0);
1343   esd->SetTimeStamp(0);
1344   esd->SetEventType(0);
1345   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1346   if (eventHeader){
1347
1348     const UInt_t *id = eventHeader->GetP("Id");
1349     esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1350     esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1351     esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1352
1353     esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1354     esd->SetEventType((eventHeader->Get("Type")));
1355   }
1356
1357   return kTRUE;
1358 }
1359
1360
1361 //_____________________________________________________________________________
1362 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1363 {
1364 // check whether detName is contained in detectors
1365 // if yes, it is removed from detectors
1366
1367   // check if all detectors are selected
1368   if ((detectors.CompareTo("ALL") == 0) ||
1369       detectors.BeginsWith("ALL ") ||
1370       detectors.EndsWith(" ALL") ||
1371       detectors.Contains(" ALL ")) {
1372     detectors = "ALL";
1373     return kTRUE;
1374   }
1375
1376   // search for the given detector
1377   Bool_t result = kFALSE;
1378   if ((detectors.CompareTo(detName) == 0) ||
1379       detectors.BeginsWith(detName+" ") ||
1380       detectors.EndsWith(" "+detName) ||
1381       detectors.Contains(" "+detName+" ")) {
1382     detectors.ReplaceAll(detName, "");
1383     result = kTRUE;
1384   }
1385
1386   // clean up the detectors string
1387   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1388   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1389   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1390
1391   return result;
1392 }
1393
1394 //_____________________________________________________________________________
1395 Bool_t AliReconstruction::InitRunLoader()
1396 {
1397 // get or create the run loader
1398
1399   if (gAlice) delete gAlice;
1400   gAlice = NULL;
1401
1402   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1403     // load all base libraries to get the loader classes
1404     TString libs = gSystem->GetLibraries();
1405     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1406       TString detName = fgkDetectorName[iDet];
1407       if (detName == "HLT") continue;
1408       if (libs.Contains("lib" + detName + "base.so")) continue;
1409       gSystem->Load("lib" + detName + "base.so");
1410     }
1411     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1412     if (!fRunLoader) {
1413       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1414       CleanUp();
1415       return kFALSE;
1416     }
1417     fRunLoader->CdGAFile();
1418     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1419       if (fRunLoader->LoadgAlice() == 0) {
1420         gAlice = fRunLoader->GetAliRun();
1421         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1422       }
1423     }
1424     if (!gAlice && !fRawReader) {
1425       AliError(Form("no gAlice object found in file %s",
1426                     fGAliceFileName.Data()));
1427       CleanUp();
1428       return kFALSE;
1429     }
1430
1431   } else {               // galice.root does not exist
1432     if (!fRawReader) {
1433       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1434       CleanUp();
1435       return kFALSE;
1436     }
1437     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1438                                     AliConfig::GetDefaultEventFolderName(),
1439                                     "recreate");
1440     if (!fRunLoader) {
1441       AliError(Form("could not create run loader in file %s",
1442                     fGAliceFileName.Data()));
1443       CleanUp();
1444       return kFALSE;
1445     }
1446     fRunLoader->MakeTree("E");
1447     Int_t iEvent = 0;
1448     while (fRawReader->NextEvent()) {
1449       fRunLoader->SetEventNumber(iEvent);
1450       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1451                                      iEvent, iEvent);
1452       fRunLoader->MakeTree("H");
1453       fRunLoader->TreeE()->Fill();
1454       iEvent++;
1455     }
1456     fRawReader->RewindEvents();
1457     if (fNumberOfEventsPerFile > 0)
1458       fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1459     else
1460       fRunLoader->SetNumberOfEventsPerFile(iEvent);
1461     fRunLoader->WriteHeader("OVERWRITE");
1462     fRunLoader->CdGAFile();
1463     fRunLoader->Write(0, TObject::kOverwrite);
1464 //    AliTracker::SetFieldMap(???);
1465   }
1466
1467   return kTRUE;
1468 }
1469
1470 //_____________________________________________________________________________
1471 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1472 {
1473 // get the reconstructor object and the loader for a detector
1474
1475   if (fReconstructor[iDet]) return fReconstructor[iDet];
1476
1477   // load the reconstructor object
1478   TPluginManager* pluginManager = gROOT->GetPluginManager();
1479   TString detName = fgkDetectorName[iDet];
1480   TString recName = "Ali" + detName + "Reconstructor";
1481   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1482
1483   if (detName == "HLT") {
1484     if (!gROOT->GetClass("AliLevel3")) {
1485       gSystem->Load("libAliHLTSrc.so");
1486       gSystem->Load("libAliHLTMisc.so");
1487       gSystem->Load("libAliHLTHough.so");
1488       gSystem->Load("libAliHLTComp.so");
1489     }
1490   }
1491
1492   AliReconstructor* reconstructor = NULL;
1493   // first check if a plugin is defined for the reconstructor
1494   TPluginHandler* pluginHandler =
1495     pluginManager->FindHandler("AliReconstructor", detName);
1496   // if not, add a plugin for it
1497   if (!pluginHandler) {
1498     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1499     TString libs = gSystem->GetLibraries();
1500     if (libs.Contains("lib" + detName + "base.so") ||
1501         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1502       pluginManager->AddHandler("AliReconstructor", detName,
1503                                 recName, detName + "rec", recName + "()");
1504     } else {
1505       pluginManager->AddHandler("AliReconstructor", detName,
1506                                 recName, detName, recName + "()");
1507     }
1508     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1509   }
1510   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1511     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1512   }
1513   if (reconstructor) {
1514     TObject* obj = fOptions.FindObject(detName.Data());
1515     if (obj) reconstructor->SetOption(obj->GetTitle());
1516     reconstructor->Init(fRunLoader);
1517     fReconstructor[iDet] = reconstructor;
1518   }
1519
1520   // get or create the loader
1521   if (detName != "HLT") {
1522     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1523     if (!fLoader[iDet]) {
1524       AliConfig::Instance()
1525         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1526                                 detName, detName);
1527       // first check if a plugin is defined for the loader
1528       pluginHandler =
1529         pluginManager->FindHandler("AliLoader", detName);
1530       // if not, add a plugin for it
1531       if (!pluginHandler) {
1532         TString loaderName = "Ali" + detName + "Loader";
1533         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1534         pluginManager->AddHandler("AliLoader", detName,
1535                                   loaderName, detName + "base",
1536                                   loaderName + "(const char*, TFolder*)");
1537         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1538       }
1539       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1540         fLoader[iDet] =
1541           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1542                                                  fRunLoader->GetEventFolder());
1543       }
1544       if (!fLoader[iDet]) {   // use default loader
1545         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1546       }
1547       if (!fLoader[iDet]) {
1548         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1549         if (fStopOnError) return NULL;
1550       } else {
1551         fRunLoader->AddLoader(fLoader[iDet]);
1552         fRunLoader->CdGAFile();
1553         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1554         fRunLoader->Write(0, TObject::kOverwrite);
1555       }
1556     }
1557   }
1558
1559   return reconstructor;
1560 }
1561
1562 //_____________________________________________________________________________
1563 Bool_t AliReconstruction::CreateVertexer()
1564 {
1565 // create the vertexer
1566
1567   fVertexer = NULL;
1568   AliReconstructor* itsReconstructor = GetReconstructor(0);
1569   if (itsReconstructor) {
1570     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1571   }
1572   if (!fVertexer) {
1573     AliWarning("couldn't create a vertexer for ITS");
1574     if (fStopOnError) return kFALSE;
1575   }
1576
1577   return kTRUE;
1578 }
1579
1580 //_____________________________________________________________________________
1581 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1582 {
1583 // create the trackers
1584
1585   TString detStr = detectors;
1586   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1587     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1588     AliReconstructor* reconstructor = GetReconstructor(iDet);
1589     if (!reconstructor) continue;
1590     TString detName = fgkDetectorName[iDet];
1591     if (detName == "HLT") {
1592       fRunHLTTracking = kTRUE;
1593       continue;
1594     }
1595     if (detName == "MUON") {
1596       fRunMuonTracking = kTRUE;
1597       continue;
1598     }
1599
1600
1601     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1602     if (!fTracker[iDet] && (iDet < 7)) {
1603       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1604       if (fStopOnError) return kFALSE;
1605     }
1606   }
1607
1608   return kTRUE;
1609 }
1610
1611 //_____________________________________________________________________________
1612 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1613 {
1614 // delete trackers and the run loader and close and delete the file
1615
1616   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1617     delete fReconstructor[iDet];
1618     fReconstructor[iDet] = NULL;
1619     fLoader[iDet] = NULL;
1620     delete fTracker[iDet];
1621     fTracker[iDet] = NULL;
1622   }
1623   delete fVertexer;
1624   fVertexer = NULL;
1625   delete fDiamondProfile;
1626   fDiamondProfile = NULL;
1627
1628   delete fRunLoader;
1629   fRunLoader = NULL;
1630   delete fRawReader;
1631   fRawReader = NULL;
1632
1633   if (file) {
1634     file->Close();
1635     delete file;
1636   }
1637
1638   if (fileOld) {
1639     fileOld->Close();
1640     delete fileOld;
1641     gSystem->Unlink("AliESDs.old.root");
1642   }
1643 }
1644
1645
1646 //_____________________________________________________________________________
1647 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1648 {
1649 // read the ESD event from a file
1650
1651   if (!esd) return kFALSE;
1652   char fileName[256];
1653   sprintf(fileName, "ESD_%d.%d_%s.root",
1654           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1655   if (gSystem->AccessPathName(fileName)) return kFALSE;
1656
1657   AliInfo(Form("reading ESD from file %s", fileName));
1658   AliDebug(1, Form("reading ESD from file %s", fileName));
1659   TFile* file = TFile::Open(fileName);
1660   if (!file || !file->IsOpen()) {
1661     AliError(Form("opening %s failed", fileName));
1662     delete file;
1663     return kFALSE;
1664   }
1665
1666   gROOT->cd();
1667   delete esd;
1668   esd = (AliESD*) file->Get("ESD");
1669   file->Close();
1670   delete file;
1671   return kTRUE;
1672 }
1673
1674 //_____________________________________________________________________________
1675 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1676 {
1677 // write the ESD event to a file
1678
1679   if (!esd) return;
1680   char fileName[256];
1681   sprintf(fileName, "ESD_%d.%d_%s.root",
1682           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1683
1684   AliDebug(1, Form("writing ESD to file %s", fileName));
1685   TFile* file = TFile::Open(fileName, "recreate");
1686   if (!file || !file->IsOpen()) {
1687     AliError(Form("opening %s failed", fileName));
1688   } else {
1689     esd->Write("ESD");
1690     file->Close();
1691   }
1692   delete file;
1693 }
1694
1695
1696
1697
1698 //_____________________________________________________________________________
1699 void AliReconstruction::CreateTag(TFile* file)
1700 {
1701   //GRP
1702   Float_t lhcLuminosity = 0.0;
1703   TString lhcState = "test";
1704   UInt_t detectorMask = 0;
1705
1706   /////////////
1707   //muon code//
1708   ////////////
1709   Double_t fMUONMASS = 0.105658369;
1710   //Variables
1711   Double_t fX,fY,fZ ;
1712   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1713   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1714   Int_t fCharge;
1715   TLorentzVector fEPvector;
1716
1717   Float_t fZVertexCut = 10.0;
1718   Float_t fRhoVertexCut = 2.0;
1719
1720   Float_t fLowPtCut = 1.0;
1721   Float_t fHighPtCut = 3.0;
1722   Float_t fVeryHighPtCut = 10.0;
1723   ////////////
1724
1725   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1726
1727   // Creates the tags for all the events in a given ESD file
1728   Int_t ntrack;
1729   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1730   Int_t nPos, nNeg, nNeutr;
1731   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1732   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1733   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1734   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1735   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1736   Int_t fVertexflag;
1737   Int_t iRunNumber = 0;
1738   TString fVertexName("default");
1739
1740   AliRunTag *tag = new AliRunTag();
1741   AliEventTag *evTag = new AliEventTag();
1742   TTree ttag("T","A Tree with event tags");
1743   TBranch * btag = ttag.Branch("AliTAG", &tag);
1744   btag->SetCompressionLevel(9);
1745
1746   AliInfo(Form("Creating the tags......."));
1747
1748   if (!file || !file->IsOpen()) {
1749     AliError(Form("opening failed"));
1750     delete file;
1751     return ;
1752   }
1753   Int_t lastEvent = 0;
1754   TTree *t = (TTree*) file->Get("esdTree");
1755   TBranch * b = t->GetBranch("ESD");
1756   AliESD *esd = 0;
1757   b->SetAddress(&esd);
1758
1759   b->GetEntry(fFirstEvent);
1760   Int_t iInitRunNumber = esd->GetRunNumber();
1761
1762   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1763   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1764   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1765     ntrack = 0;
1766     nPos = 0;
1767     nNeg = 0;
1768     nNeutr =0;
1769     nK0s = 0;
1770     nNeutrons = 0;
1771     nPi0s = 0;
1772     nGamas = 0;
1773     nProtons = 0;
1774     nKaons = 0;
1775     nPions = 0;
1776     nMuons = 0;
1777     nElectrons = 0;
1778     nCh1GeV = 0;
1779     nCh3GeV = 0;
1780     nCh10GeV = 0;
1781     nMu1GeV = 0;
1782     nMu3GeV = 0;
1783     nMu10GeV = 0;
1784     nEl1GeV = 0;
1785     nEl3GeV = 0;
1786     nEl10GeV = 0;
1787     maxPt = .0;
1788     meanPt = .0;
1789     totalP = .0;
1790     fVertexflag = 0;
1791
1792     b->GetEntry(iEventNumber);
1793     iRunNumber = esd->GetRunNumber();
1794     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1795     const AliESDVertex * vertexIn = esd->GetVertex();
1796     if (!vertexIn) AliError("ESD has not defined vertex.");
1797     if (vertexIn) fVertexName = vertexIn->GetName();
1798     if(fVertexName != "default") fVertexflag = 1;
1799     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1800       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1801       UInt_t status = esdTrack->GetStatus();
1802
1803       //select only tracks with ITS refit
1804       if ((status&AliESDtrack::kITSrefit)==0) continue;
1805       //select only tracks with TPC refit
1806       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1807
1808       //select only tracks with the "combined PID"
1809       if ((status&AliESDtrack::kESDpid)==0) continue;
1810       Double_t p[3];
1811       esdTrack->GetPxPyPz(p);
1812       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1813       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1814       totalP += momentum;
1815       meanPt += fPt;
1816       if(fPt > maxPt) maxPt = fPt;
1817
1818       if(esdTrack->GetSign() > 0) {
1819         nPos++;
1820         if(fPt > fLowPtCut) nCh1GeV++;
1821         if(fPt > fHighPtCut) nCh3GeV++;
1822         if(fPt > fVeryHighPtCut) nCh10GeV++;
1823       }
1824       if(esdTrack->GetSign() < 0) {
1825         nNeg++;
1826         if(fPt > fLowPtCut) nCh1GeV++;
1827         if(fPt > fHighPtCut) nCh3GeV++;
1828         if(fPt > fVeryHighPtCut) nCh10GeV++;
1829       }
1830       if(esdTrack->GetSign() == 0) nNeutr++;
1831
1832       //PID
1833       Double_t prob[5];
1834       esdTrack->GetESDpid(prob);
1835
1836       Double_t rcc = 0.0;
1837       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1838       if(rcc == 0.0) continue;
1839       //Bayes' formula
1840       Double_t w[5];
1841       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1842
1843       //protons
1844       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1845       //kaons
1846       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1847       //pions
1848       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1849       //electrons
1850       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1851         nElectrons++;
1852         if(fPt > fLowPtCut) nEl1GeV++;
1853         if(fPt > fHighPtCut) nEl3GeV++;
1854         if(fPt > fVeryHighPtCut) nEl10GeV++;
1855       }
1856       ntrack++;
1857     }//track loop
1858
1859     /////////////
1860     //muon code//
1861     ////////////
1862     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1863     // loop over all reconstructed tracks (also first track of combination)
1864     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1865       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1866       if (muonTrack == 0x0) continue;
1867
1868       // Coordinates at vertex
1869       fZ = muonTrack->GetZ();
1870       fY = muonTrack->GetBendingCoor();
1871       fX = muonTrack->GetNonBendingCoor();
1872
1873       fThetaX = muonTrack->GetThetaX();
1874       fThetaY = muonTrack->GetThetaY();
1875
1876       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1877       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1878       fPxRec = fPzRec * TMath::Tan(fThetaX);
1879       fPyRec = fPzRec * TMath::Tan(fThetaY);
1880       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1881
1882       //ChiSquare of the track if needed
1883       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1884       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1885       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1886
1887       // total number of muons inside a vertex cut
1888       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1889         nMuons++;
1890         if(fEPvector.Pt() > fLowPtCut) {
1891           nMu1GeV++;
1892           if(fEPvector.Pt() > fHighPtCut) {
1893             nMu3GeV++;
1894             if (fEPvector.Pt() > fVeryHighPtCut) {
1895               nMu10GeV++;
1896             }
1897           }
1898         }
1899       }
1900     }//muon track loop
1901
1902     // Fill the event tags
1903     if(ntrack != 0)
1904       meanPt = meanPt/ntrack;
1905
1906     evTag->SetEventId(iEventNumber+1);
1907     if (vertexIn) {
1908       evTag->SetVertexX(vertexIn->GetXv());
1909       evTag->SetVertexY(vertexIn->GetYv());
1910       evTag->SetVertexZ(vertexIn->GetZv());
1911       evTag->SetVertexZError(vertexIn->GetZRes());
1912     }
1913     evTag->SetVertexFlag(fVertexflag);
1914
1915     evTag->SetT0VertexZ(esd->GetT0zVertex());
1916
1917     evTag->SetTriggerMask(esd->GetTriggerMask());
1918     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1919
1920     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1921     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1922     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1923     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1924     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1925     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1926
1927
1928     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1929     evTag->SetNumOfPosTracks(nPos);
1930     evTag->SetNumOfNegTracks(nNeg);
1931     evTag->SetNumOfNeutrTracks(nNeutr);
1932
1933     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1934     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1935     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1936     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1937
1938     evTag->SetNumOfProtons(nProtons);
1939     evTag->SetNumOfKaons(nKaons);
1940     evTag->SetNumOfPions(nPions);
1941     evTag->SetNumOfMuons(nMuons);
1942     evTag->SetNumOfElectrons(nElectrons);
1943     evTag->SetNumOfPhotons(nGamas);
1944     evTag->SetNumOfPi0s(nPi0s);
1945     evTag->SetNumOfNeutrons(nNeutrons);
1946     evTag->SetNumOfKaon0s(nK0s);
1947
1948     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1949     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1950     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1951     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1952     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1953     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1954     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1955     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1956     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1957
1958     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1959     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1960
1961     evTag->SetTotalMomentum(totalP);
1962     evTag->SetMeanPt(meanPt);
1963     evTag->SetMaxPt(maxPt);
1964
1965     tag->SetLHCTag(lhcLuminosity,lhcState);
1966     tag->SetDetectorTag(detectorMask);
1967
1968     tag->SetRunId(iInitRunNumber);
1969     tag->AddEventTag(*evTag);
1970   }
1971   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
1972   else lastEvent = fLastEvent;
1973
1974   ttag.Fill();
1975   tag->Clear();
1976
1977   char fileName[256];
1978   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1979           tag->GetRunId(),fFirstEvent,lastEvent );
1980   AliInfo(Form("writing tags to file %s", fileName));
1981   AliDebug(1, Form("writing tags to file %s", fileName));
1982
1983   TFile* ftag = TFile::Open(fileName, "recreate");
1984   ftag->cd();
1985   ttag.Write();
1986   ftag->Close();
1987   file->cd();
1988   delete tag;
1989   delete evTag;
1990 }
1991
1992 //_____________________________________________________________________________
1993 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
1994 {
1995   // write all files from the given esd file to an aod file
1996
1997   // create an AliAOD object
1998   AliAODEvent *aod = new AliAODEvent();
1999   aod->CreateStdContent();
2000
2001   // go to the file
2002   aodFile->cd();
2003
2004   // create the tree
2005   TTree *aodTree = new TTree("AOD", "AliAOD tree");
2006   aodTree->Branch(aod->GetList());
2007
2008   // connect to ESD
2009   TTree *t = (TTree*) esdFile->Get("esdTree");
2010   TBranch *b = t->GetBranch("ESD");
2011   AliESD *esd = 0;
2012   b->SetAddress(&esd);
2013
2014   Int_t nEvents = b->GetEntries();
2015
2016   // set arrays and pointers
2017   Float_t posF[3];
2018   Double_t pos[3];
2019   Double_t p[3];
2020   Double_t covVtx[6];
2021   Double_t covTr[21];
2022   Double_t pid[10];
2023
2024   // loop over events and fill them
2025   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2026     b->GetEntry(iEvent);
2027
2028     // Multiplicity information needed by the header (to be revised!)
2029     Int_t nTracks   = esd->GetNumberOfTracks();
2030     Int_t nPosTracks = 0;
2031     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2032       if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2033
2034     // create the header
2035     aod->AddHeader(new AliAODHeader(esd->GetRunNumber(),
2036                                     esd->GetBunchCrossNumber(),
2037                                     esd->GetOrbitNumber(),
2038                                     esd->GetPeriodNumber(),
2039                                     nTracks,
2040                                     nPosTracks,
2041                                     nTracks-nPosTracks,
2042                                     esd->GetMagneticField(),
2043                                     -999., // fill muon magnetic field
2044                                     -999., // centrality; to be filled, still
2045                                     esd->GetZDCN1Energy(),
2046                                     esd->GetZDCP1Energy(),
2047                                     esd->GetZDCN2Energy(),
2048                                     esd->GetZDCP2Energy(),
2049                                     esd->GetZDCEMEnergy(),
2050                                     esd->GetTriggerMask(),
2051                                     esd->GetTriggerCluster(),
2052                                     esd->GetEventType()));
2053
2054     Int_t nV0s      = esd->GetNumberOfV0s();
2055     Int_t nCascades = esd->GetNumberOfCascades();
2056     Int_t nKinks    = esd->GetNumberOfKinks();
2057     Int_t nVertices = nV0s + nCascades + nKinks;
2058
2059     aod->ResetStd(nTracks, nVertices);
2060     AliAODTrack *aodTrack;
2061
2062
2063     // Array to take into account the tracks already added to the AOD
2064     Bool_t * usedTrack = NULL;
2065     if (nTracks>0) {
2066       usedTrack = new Bool_t[nTracks];
2067       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2068     }
2069     // Array to take into account the V0s already added to the AOD
2070     Bool_t * usedV0 = NULL;
2071     if (nV0s>0) {
2072       usedV0 = new Bool_t[nV0s];
2073       for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2074     }
2075     // Array to take into account the kinks already added to the AOD
2076     Bool_t * usedKink = NULL;
2077     if (nKinks>0) {
2078       usedKink = new Bool_t[nKinks];
2079       for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2080     }
2081
2082     // Access to the AOD container of vertices
2083     TClonesArray &vertices = *(aod->GetVertices());
2084     Int_t jVertices=0;
2085
2086     // Access to the AOD container of tracks
2087     TClonesArray &tracks = *(aod->GetTracks());
2088     Int_t jTracks=0;
2089
2090     // Add primary vertex. The primary tracks will be defined
2091     // after the loops on the composite objects (V0, cascades, kinks)
2092     const AliESDVertex *vtx = esd->GetPrimaryVertex();
2093
2094     vtx->GetXYZ(pos); // position
2095     vtx->GetCovMatrix(covVtx); //covariance matrix
2096
2097     AliAODVertex * primary = new(vertices[jVertices++])
2098       AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2099
2100     // Create vertices starting from the most complex objects
2101
2102     // Cascades
2103     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2104       AliESDcascade *cascade = esd->GetCascade(nCascade);
2105
2106       cascade->GetXYZ(pos[0], pos[1], pos[2]);
2107       cascade->GetPosCovXi(covVtx);
2108
2109       // Add the cascade vertex
2110       AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2111                                                                         covVtx,
2112                                                                         cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2113                                                                         primary,
2114                                                                         AliAODVertex::kCascade);
2115
2116       primary->AddDaughter(vcascade);
2117
2118       // Add the V0 from the cascade. The ESD class have to be optimized...
2119       // Now we have to search for the corresponding Vo in the list of V0s
2120       // using the indeces of the positive and negative tracks
2121
2122       Int_t posFromV0 = cascade->GetPindex();
2123       Int_t negFromV0 = cascade->GetNindex();
2124
2125
2126       AliESDv0 * v0 = 0x0;
2127       Int_t indV0 = -1;
2128
2129       for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2130
2131         v0 = esd->GetV0(iV0);
2132         Int_t posV0 = v0->GetPindex();
2133         Int_t negV0 = v0->GetNindex();
2134
2135         if (posV0==posFromV0 && negV0==negFromV0) {
2136           indV0 = iV0;
2137           break;
2138         }
2139       }
2140
2141       AliAODVertex * vV0FromCascade = 0x0;
2142
2143       if (indV0>-1 && !usedV0[indV0] ) {
2144
2145         // the V0 exists in the array of V0s and is not used
2146
2147         usedV0[indV0] = kTRUE;
2148
2149         v0->GetXYZ(pos[0], pos[1], pos[2]);
2150         v0->GetPosCov(covVtx);
2151
2152         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2153                                                                  covVtx,
2154                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2155                                                                  vcascade,
2156                                                                  AliAODVertex::kV0);
2157       } else {
2158
2159         // the V0 doesn't exist in the array of V0s or was used
2160         cerr << "Error: event " << iEvent << " cascade " << nCascade
2161              << " The V0 " << indV0
2162              << " doesn't exist in the array of V0s or was used!" << endl;
2163
2164         cascade->GetXYZ(pos[0], pos[1], pos[2]);
2165         cascade->GetPosCov(covVtx);
2166
2167         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2168                                                                  covVtx,
2169                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2170                                                                  vcascade,
2171                                                                  AliAODVertex::kV0);
2172         vcascade->AddDaughter(vV0FromCascade);
2173       }
2174
2175       // Add the positive tracks from the V0
2176
2177       if (! usedTrack[posFromV0]) {
2178
2179         usedTrack[posFromV0] = kTRUE;
2180
2181         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2182         esdTrack->GetPxPyPz(p);
2183         esdTrack->GetXYZ(pos);
2184         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2185         esdTrack->GetESDpid(pid);
2186
2187         vV0FromCascade->AddDaughter(aodTrack =
2188                                     new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2189                                            esdTrack->GetLabel(),
2190                                            p,
2191                                            kTRUE,
2192                                            pos,
2193                                            kFALSE,
2194                                            covTr,
2195                                            (Short_t)esdTrack->GetSign(),
2196                                            esdTrack->GetITSClusterMap(),
2197                                            pid,
2198                                            vV0FromCascade,
2199                                            kTRUE,  // check if this is right
2200                                            kFALSE, // check if this is right
2201                                            AliAODTrack::kSecondary)
2202                 );
2203         aodTrack->ConvertAliPIDtoAODPID();
2204       }
2205       else {
2206         cerr << "Error: event " << iEvent << " cascade " << nCascade
2207              << " track " << posFromV0 << " has already been used!" << endl;
2208       }
2209
2210       // Add the negative tracks from the V0
2211
2212       if (!usedTrack[negFromV0]) {
2213
2214         usedTrack[negFromV0] = kTRUE;
2215
2216         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2217         esdTrack->GetPxPyPz(p);
2218         esdTrack->GetXYZ(pos);
2219         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2220         esdTrack->GetESDpid(pid);
2221
2222         vV0FromCascade->AddDaughter(aodTrack =
2223                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2224                                            esdTrack->GetLabel(),
2225                                            p,
2226                                            kTRUE,
2227                                            pos,
2228                                            kFALSE,
2229                                            covTr,
2230                                            (Short_t)esdTrack->GetSign(),
2231                                            esdTrack->GetITSClusterMap(),
2232                                            pid,
2233                                            vV0FromCascade,
2234                                            kTRUE,  // check if this is right
2235                                            kFALSE, // check if this is right
2236                                            AliAODTrack::kSecondary)
2237                 );
2238         aodTrack->ConvertAliPIDtoAODPID();
2239       }
2240       else {
2241         cerr << "Error: event " << iEvent << " cascade " << nCascade
2242              << " track " << negFromV0 << " has already been used!" << endl;
2243       }
2244
2245       // Add the bachelor track from the cascade
2246
2247       Int_t bachelor = cascade->GetBindex();
2248
2249       if(!usedTrack[bachelor]) {
2250
2251         usedTrack[bachelor] = kTRUE;
2252
2253         AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2254         esdTrack->GetPxPyPz(p);
2255         esdTrack->GetXYZ(pos);
2256         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2257         esdTrack->GetESDpid(pid);
2258
2259         vcascade->AddDaughter(aodTrack =
2260                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2261                                            esdTrack->GetLabel(),
2262                                            p,
2263                                            kTRUE,
2264                                            pos,
2265                                            kFALSE,
2266                                            covTr,
2267                                            (Short_t)esdTrack->GetSign(),
2268                                            esdTrack->GetITSClusterMap(),
2269                                            pid,
2270                                            vcascade,
2271                                            kTRUE,  // check if this is right
2272                                            kFALSE, // check if this is right
2273                                            AliAODTrack::kSecondary)
2274                 );
2275         aodTrack->ConvertAliPIDtoAODPID();
2276      }
2277       else {
2278         cerr << "Error: event " << iEvent << " cascade " << nCascade
2279              << " track " << bachelor << " has already been used!" << endl;
2280       }
2281
2282       // Add the primary track of the cascade (if any)
2283
2284     } // end of the loop on cascades
2285
2286     // V0s
2287
2288     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2289
2290       if (usedV0[nV0]) continue; // skip if aready added to the AOD
2291
2292       AliESDv0 *v0 = esd->GetV0(nV0);
2293
2294       v0->GetXYZ(pos[0], pos[1], pos[2]);
2295       v0->GetPosCov(covVtx);
2296
2297       AliAODVertex * vV0 =
2298         new(vertices[jVertices++]) AliAODVertex(pos,
2299                                                 covVtx,
2300                                                 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2301                                                 primary,
2302                                                 AliAODVertex::kV0);
2303       primary->AddDaughter(vV0);
2304
2305       Int_t posFromV0 = v0->GetPindex();
2306       Int_t negFromV0 = v0->GetNindex();
2307
2308       // Add the positive tracks from the V0
2309
2310       if (!usedTrack[posFromV0]) {
2311
2312         usedTrack[posFromV0] = kTRUE;
2313
2314         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2315         esdTrack->GetPxPyPz(p);
2316         esdTrack->GetXYZ(pos);
2317         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2318         esdTrack->GetESDpid(pid);
2319
2320         vV0->AddDaughter(aodTrack =
2321                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2322                                            esdTrack->GetLabel(),
2323                                            p,
2324                                            kTRUE,
2325                                            pos,
2326                                            kFALSE,
2327                                            covTr,
2328                                            (Short_t)esdTrack->GetSign(),
2329                                            esdTrack->GetITSClusterMap(),
2330                                            pid,
2331                                            vV0,
2332                                            kTRUE,  // check if this is right
2333                                            kFALSE, // check if this is right
2334                                            AliAODTrack::kSecondary)
2335                 );
2336         aodTrack->ConvertAliPIDtoAODPID();
2337       }
2338       else {
2339         cerr << "Error: event " << iEvent << " V0 " << nV0
2340              << " track " << posFromV0 << " has already been used!" << endl;
2341       }
2342
2343       // Add the negative tracks from the V0
2344
2345       if (!usedTrack[negFromV0]) {
2346
2347         usedTrack[negFromV0] = kTRUE;
2348
2349         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2350         esdTrack->GetPxPyPz(p);
2351         esdTrack->GetXYZ(pos);
2352         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2353         esdTrack->GetESDpid(pid);
2354
2355         vV0->AddDaughter(aodTrack =
2356                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2357                                            esdTrack->GetLabel(),
2358                                            p,
2359                                            kTRUE,
2360                                            pos,
2361                                            kFALSE,
2362                                            covTr,
2363                                            (Short_t)esdTrack->GetSign(),
2364                                            esdTrack->GetITSClusterMap(),
2365                                            pid,
2366                                            vV0,
2367                                            kTRUE,  // check if this is right
2368                                            kFALSE, // check if this is right
2369                                            AliAODTrack::kSecondary)
2370                 );
2371         aodTrack->ConvertAliPIDtoAODPID();
2372       }
2373       else {
2374         cerr << "Error: event " << iEvent << " V0 " << nV0
2375              << " track " << negFromV0 << " has already been used!" << endl;
2376       }
2377
2378     } // end of the loop on V0s
2379
2380     // Kinks: it is a big mess the access to the information in the kinks
2381     // The loop is on the tracks in order to find the mother and daugther of each kink
2382
2383
2384     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2385
2386
2387       AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2388
2389       Int_t ikink = esdTrack->GetKinkIndex(0);
2390
2391       if (ikink) {
2392         // Negative kink index: mother, positive: daughter
2393
2394         // Search for the second track of the kink
2395
2396         for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2397
2398           AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2399
2400           Int_t jkink = esdTrack1->GetKinkIndex(0);
2401
2402           if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2403
2404             // The two tracks are from the same kink
2405
2406             if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2407
2408             Int_t imother = -1;
2409             Int_t idaughter = -1;
2410
2411             if (ikink<0 && jkink>0) {
2412
2413               imother = iTrack;
2414               idaughter = jTrack;
2415             }
2416             else if (ikink>0 && jkink<0) {
2417
2418               imother = jTrack;
2419               idaughter = iTrack;
2420             }
2421             else {
2422               cerr << "Error: Wrong combination of kink indexes: "
2423               << ikink << " " << jkink << endl;
2424               continue;
2425             }
2426
2427             // Add the mother track
2428
2429             AliAODTrack * mother = NULL;
2430
2431             if (!usedTrack[imother]) {
2432
2433               usedTrack[imother] = kTRUE;
2434
2435               AliESDtrack *esdTrack = esd->GetTrack(imother);
2436               esdTrack->GetPxPyPz(p);
2437               esdTrack->GetXYZ(pos);
2438               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2439               esdTrack->GetESDpid(pid);
2440
2441               mother =
2442                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2443                                            esdTrack->GetLabel(),
2444                                            p,
2445                                            kTRUE,
2446                                            pos,
2447                                            kFALSE,
2448                                            covTr,
2449                                            (Short_t)esdTrack->GetSign(),
2450                                            esdTrack->GetITSClusterMap(),
2451                                            pid,
2452                                            primary,
2453                                            kTRUE, // check if this is right
2454                                            kTRUE, // check if this is right
2455                                            AliAODTrack::kPrimary);
2456               primary->AddDaughter(mother);
2457               mother->ConvertAliPIDtoAODPID();
2458             }
2459             else {
2460               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2461               << " track " << imother << " has already been used!" << endl;
2462             }
2463
2464             // Add the kink vertex
2465             AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2466
2467             AliAODVertex * vkink =
2468             new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2469                                                     NULL,
2470                                                     0.,
2471                                                     mother,
2472                                                     AliAODVertex::kKink);
2473             // Add the daughter track
2474
2475             AliAODTrack * daughter = NULL;
2476
2477             if (!usedTrack[idaughter]) {
2478
2479               usedTrack[idaughter] = kTRUE;
2480
2481               AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2482               esdTrack->GetPxPyPz(p);
2483               esdTrack->GetXYZ(pos);
2484               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2485               esdTrack->GetESDpid(pid);
2486
2487               daughter =
2488                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2489                                            esdTrack->GetLabel(),
2490                                            p,
2491                                            kTRUE,
2492                                            pos,
2493                                            kFALSE,
2494                                            covTr,
2495                                            (Short_t)esdTrack->GetSign(),
2496                                            esdTrack->GetITSClusterMap(),
2497                                            pid,
2498                                            vkink,
2499                                            kTRUE, // check if this is right
2500                                            kTRUE, // check if this is right
2501                                            AliAODTrack::kPrimary);
2502               vkink->AddDaughter(daughter);
2503               daughter->ConvertAliPIDtoAODPID();
2504             }
2505             else {
2506               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2507               << " track " << idaughter << " has already been used!" << endl;
2508             }
2509
2510
2511           }
2512         }
2513
2514       }
2515
2516     }
2517
2518
2519     // Tracks (primary and orphan)
2520
2521     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2522
2523
2524       if (usedTrack[nTrack]) continue;
2525
2526       AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2527       esdTrack->GetPxPyPz(p);
2528       esdTrack->GetXYZ(pos);
2529       esdTrack->GetCovarianceXYZPxPyPz(covTr);
2530       esdTrack->GetESDpid(pid);
2531
2532       Float_t impactXY, impactZ;
2533
2534       esdTrack->GetImpactParameters(impactXY,impactZ);
2535
2536       if (impactXY<3) {
2537         // track inside the beam pipe
2538
2539         primary->AddDaughter(aodTrack =
2540             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2541                                          esdTrack->GetLabel(),
2542                                          p,
2543                                          kTRUE,
2544                                          pos,
2545                                          kFALSE,
2546                                          covTr,
2547                                          (Short_t)esdTrack->GetSign(),
2548                                          esdTrack->GetITSClusterMap(),
2549                                          pid,
2550                                          primary,
2551                                          kTRUE, // check if this is right
2552                                          kTRUE, // check if this is right
2553                                          AliAODTrack::kPrimary)
2554             );
2555         aodTrack->ConvertAliPIDtoAODPID();
2556       }
2557       else {
2558         // outside the beam pipe: orphan track
2559             aodTrack =
2560             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2561                                          esdTrack->GetLabel(),
2562                                          p,
2563                                          kTRUE,
2564                                          pos,
2565                                          kFALSE,
2566                                          covTr,
2567                                          (Short_t)esdTrack->GetSign(),
2568                                          esdTrack->GetITSClusterMap(),
2569                                          pid,
2570                                          NULL,
2571                                          kFALSE, // check if this is right
2572                                          kFALSE, // check if this is right
2573                                          AliAODTrack::kOrphan);
2574             aodTrack->ConvertAliPIDtoAODPID();
2575       }
2576     } // end of loop on tracks
2577
2578     // muon tracks
2579     Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2580     for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2581
2582       AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2583       p[0] = esdMuTrack->Px();
2584       p[1] = esdMuTrack->Py();
2585       p[2] = esdMuTrack->Pz();
2586       pos[0] = primary->GetX();
2587       pos[1] = primary->GetY();
2588       pos[2] = primary->GetZ();
2589
2590       // has to be changed once the muon pid is provided by the ESD
2591       for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2592
2593       primary->AddDaughter(
2594           new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2595                                              0, // no label provided
2596                                              p,
2597                                              kTRUE,
2598                                              pos,
2599                                              kFALSE,
2600                                              NULL, // no covariance matrix provided
2601                                              (Short_t)-99, // no charge provided
2602                                              0, // no ITSClusterMap
2603                                              pid,
2604                                              primary,
2605                                              kTRUE,  // check if this is right
2606                                              kTRUE,  // not used for vertex fit
2607                                              AliAODTrack::kPrimary)
2608           );
2609     }
2610
2611     // Access to the AOD container of clusters
2612     TClonesArray &clusters = *(aod->GetClusters());
2613     Int_t jClusters=0;
2614
2615     // Calo Clusters
2616     Int_t nClusters    = esd->GetNumberOfCaloClusters();
2617
2618     for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2619
2620       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2621
2622       Int_t id = cluster->GetID();
2623       Int_t label = -1;
2624       Float_t energy = cluster->GetClusterEnergy();
2625       cluster->GetGlobalPosition(posF);
2626       AliAODVertex *prodVertex = primary;
2627       AliAODTrack *primTrack = NULL;
2628       Char_t ttype=AliAODCluster::kUndef;
2629
2630       if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2631       else if (cluster->IsEMCAL()) {
2632
2633         if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2634           ttype = AliAODCluster::kEMCALPseudoCluster;
2635         else
2636           ttype = AliAODCluster::kEMCALClusterv1;
2637
2638       }
2639
2640       new(clusters[jClusters++]) AliAODCluster(id,
2641                                                label,
2642                                                energy,
2643                                                pos,
2644                                                NULL, // no covariance matrix provided
2645                                                NULL, // no pid for clusters provided
2646                                                prodVertex,
2647                                                primTrack,
2648                                                ttype);
2649
2650     } // end of loop on calo clusters
2651
2652     delete [] usedTrack;
2653     delete [] usedV0;
2654     delete [] usedKink;
2655
2656     // fill the tree for this event
2657     aodTree->Fill();
2658   } // end of event loop
2659
2660   aodTree->GetUserInfo()->Add(aod);
2661
2662   // close ESD file
2663   esdFile->Close();
2664
2665   // write the tree to the specified file
2666   aodFile = aodTree->GetCurrentFile();
2667   aodFile->cd();
2668   aodTree->Write();
2669
2670   return;
2671 }
2672
2673
2674 void AliReconstruction::WriteAlignmentData(AliESD* esd)
2675 {
2676   // Write space-points which are then used in the alignment procedures
2677   // For the moment only ITS, TRD and TPC
2678
2679   // Load TOF clusters
2680   if (fTracker[3]){
2681     fLoader[3]->LoadRecPoints("read");
2682     TTree* tree = fLoader[3]->TreeR();
2683     if (!tree) {
2684       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2685       return;
2686     }
2687     fTracker[3]->LoadClusters(tree);
2688   }
2689   Int_t ntracks = esd->GetNumberOfTracks();
2690   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2691     {
2692       AliESDtrack *track = esd->GetTrack(itrack);
2693       Int_t nsp = 0;
2694       Int_t idx[200];
2695       for (Int_t iDet = 3; iDet >= 0; iDet--)
2696         nsp += track->GetNcls(iDet);
2697       if (nsp) {
2698         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2699         track->SetTrackPointArray(sp);
2700         Int_t isptrack = 0;
2701         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2702           AliTracker *tracker = fTracker[iDet];
2703           if (!tracker) continue;
2704           Int_t nspdet = track->GetNcls(iDet);
2705           if (nspdet <= 0) continue;
2706           track->GetClusters(iDet,idx);
2707           AliTrackPoint p;
2708           Int_t isp = 0;
2709           Int_t isp2 = 0;
2710           while (isp < nspdet) {
2711             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2712             const Int_t kNTPCmax = 159;
2713             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2714             if (!isvalid) continue;
2715             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2716           }
2717         }
2718       }
2719     }
2720   if (fTracker[3]){
2721     fTracker[3]->UnloadClusters();
2722     fLoader[3]->UnloadRecPoints();
2723   }
2724 }
2725
2726 //_____________________________________________________________________________
2727 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2728 {
2729   // The method reads the raw-data error log
2730   // accumulated within the rawReader.
2731   // It extracts the raw-data errors related to
2732   // the current event and stores them into
2733   // a TClonesArray inside the esd object.
2734
2735   if (!fRawReader) return;
2736
2737   for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2738
2739     AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2740     if (!log) continue;
2741     if (iEvent != log->GetEventNumber()) continue;
2742
2743     esd->AddRawDataErrorLog(log);
2744   }
2745
2746 }