STEER/AliReconstruction.cxx

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