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   //GRP
1780   Float_t lhcLuminosity = 0.0;
1781   TString lhcState = "test";
1782   UInt_t detectorMask = 0;
1783
1784   /////////////
1785   //muon code//
1786   ////////////
1787   Double_t fMUONMASS = 0.105658369;
1788   //Variables
1789   Double_t fX,fY,fZ ;
1790   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1791   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1792   Int_t fCharge;
1793   TLorentzVector fEPvector;
1794
1795   Float_t fZVertexCut = 10.0;
1796   Float_t fRhoVertexCut = 2.0;
1797
1798   Float_t fLowPtCut = 1.0;
1799   Float_t fHighPtCut = 3.0;
1800   Float_t fVeryHighPtCut = 10.0;
1801   ////////////
1802
1803   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1804
1805   // Creates the tags for all the events in a given ESD file
1806   Int_t ntrack;
1807   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1808   Int_t nPos, nNeg, nNeutr;
1809   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1810   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1811   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1812   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1813   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1814   Int_t fVertexflag;
1815   Int_t iRunNumber = 0;
1816   TString fVertexName("default");
1817
1818   AliRunTag *tag = new AliRunTag();
1819   AliEventTag *evTag = new AliEventTag();
1820   TTree ttag("T","A Tree with event tags");
1821   TBranch * btag = ttag.Branch("AliTAG", &tag);
1822   btag->SetCompressionLevel(9);
1823
1824   AliInfo(Form("Creating the tags......."));
1825
1826   if (!file || !file->IsOpen()) {
1827     AliError(Form("opening failed"));
1828     delete file;
1829     return ;
1830   }
1831   Int_t lastEvent = 0;
1832   TTree *t = (TTree*) file->Get("esdTree");
1833   TBranch * b = t->GetBranch("ESD");
1834   AliESD *esd = 0;
1835   b->SetAddress(&esd);
1836
1837   b->GetEntry(fFirstEvent);
1838   Int_t iInitRunNumber = esd->GetRunNumber();
1839
1840   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1841   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1842   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1843     ntrack = 0;
1844     nPos = 0;
1845     nNeg = 0;
1846     nNeutr =0;
1847     nK0s = 0;
1848     nNeutrons = 0;
1849     nPi0s = 0;
1850     nGamas = 0;
1851     nProtons = 0;
1852     nKaons = 0;
1853     nPions = 0;
1854     nMuons = 0;
1855     nElectrons = 0;
1856     nCh1GeV = 0;
1857     nCh3GeV = 0;
1858     nCh10GeV = 0;
1859     nMu1GeV = 0;
1860     nMu3GeV = 0;
1861     nMu10GeV = 0;
1862     nEl1GeV = 0;
1863     nEl3GeV = 0;
1864     nEl10GeV = 0;
1865     maxPt = .0;
1866     meanPt = .0;
1867     totalP = .0;
1868     fVertexflag = 0;
1869
1870     b->GetEntry(iEventNumber);
1871     iRunNumber = esd->GetRunNumber();
1872     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1873     const AliESDVertex * vertexIn = esd->GetVertex();
1874     if (!vertexIn) AliError("ESD has not defined vertex.");
1875     if (vertexIn) fVertexName = vertexIn->GetName();
1876     if(fVertexName != "default") fVertexflag = 1;
1877     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1878       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1879       UInt_t status = esdTrack->GetStatus();
1880
1881       //select only tracks with ITS refit
1882       if ((status&AliESDtrack::kITSrefit)==0) continue;
1883       //select only tracks with TPC refit
1884       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1885
1886       //select only tracks with the "combined PID"
1887       if ((status&AliESDtrack::kESDpid)==0) continue;
1888       Double_t p[3];
1889       esdTrack->GetPxPyPz(p);
1890       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1891       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1892       totalP += momentum;
1893       meanPt += fPt;
1894       if(fPt > maxPt) maxPt = fPt;
1895
1896       if(esdTrack->GetSign() > 0) {
1897         nPos++;
1898         if(fPt > fLowPtCut) nCh1GeV++;
1899         if(fPt > fHighPtCut) nCh3GeV++;
1900         if(fPt > fVeryHighPtCut) nCh10GeV++;
1901       }
1902       if(esdTrack->GetSign() < 0) {
1903         nNeg++;
1904         if(fPt > fLowPtCut) nCh1GeV++;
1905         if(fPt > fHighPtCut) nCh3GeV++;
1906         if(fPt > fVeryHighPtCut) nCh10GeV++;
1907       }
1908       if(esdTrack->GetSign() == 0) nNeutr++;
1909
1910       //PID
1911       Double_t prob[5];
1912       esdTrack->GetESDpid(prob);
1913
1914       Double_t rcc = 0.0;
1915       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1916       if(rcc == 0.0) continue;
1917       //Bayes' formula
1918       Double_t w[5];
1919       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1920
1921       //protons
1922       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1923       //kaons
1924       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1925       //pions
1926       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1927       //electrons
1928       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1929         nElectrons++;
1930         if(fPt > fLowPtCut) nEl1GeV++;
1931         if(fPt > fHighPtCut) nEl3GeV++;
1932         if(fPt > fVeryHighPtCut) nEl10GeV++;
1933       }
1934       ntrack++;
1935     }//track loop
1936
1937     /////////////
1938     //muon code//
1939     ////////////
1940     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1941     // loop over all reconstructed tracks (also first track of combination)
1942     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1943       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1944       if (muonTrack == 0x0) continue;
1945
1946       // Coordinates at vertex
1947       fZ = muonTrack->GetZ();
1948       fY = muonTrack->GetBendingCoor();
1949       fX = muonTrack->GetNonBendingCoor();
1950
1951       fThetaX = muonTrack->GetThetaX();
1952       fThetaY = muonTrack->GetThetaY();
1953
1954       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1955       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1956       fPxRec = fPzRec * TMath::Tan(fThetaX);
1957       fPyRec = fPzRec * TMath::Tan(fThetaY);
1958       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1959
1960       //ChiSquare of the track if needed
1961       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1962       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1963       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1964
1965       // total number of muons inside a vertex cut
1966       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1967         nMuons++;
1968         if(fEPvector.Pt() > fLowPtCut) {
1969           nMu1GeV++;
1970           if(fEPvector.Pt() > fHighPtCut) {
1971             nMu3GeV++;
1972             if (fEPvector.Pt() > fVeryHighPtCut) {
1973               nMu10GeV++;
1974             }
1975           }
1976         }
1977       }
1978     }//muon track loop
1979
1980     // Fill the event tags
1981     if(ntrack != 0)
1982       meanPt = meanPt/ntrack;
1983
1984     evTag->SetEventId(iEventNumber+1);
1985     if (vertexIn) {
1986       evTag->SetVertexX(vertexIn->GetXv());
1987       evTag->SetVertexY(vertexIn->GetYv());
1988       evTag->SetVertexZ(vertexIn->GetZv());
1989       evTag->SetVertexZError(vertexIn->GetZRes());
1990     }
1991     evTag->SetVertexFlag(fVertexflag);
1992
1993     evTag->SetT0VertexZ(esd->GetT0zVertex());
1994
1995     evTag->SetTriggerMask(esd->GetTriggerMask());
1996     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1997
1998     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1999     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
2000     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
2001     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2002     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
2003     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
2004
2005
2006     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
2007     evTag->SetNumOfPosTracks(nPos);
2008     evTag->SetNumOfNegTracks(nNeg);
2009     evTag->SetNumOfNeutrTracks(nNeutr);
2010
2011     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
2012     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
2013     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
2014     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
2015
2016     evTag->SetNumOfProtons(nProtons);
2017     evTag->SetNumOfKaons(nKaons);
2018     evTag->SetNumOfPions(nPions);
2019     evTag->SetNumOfMuons(nMuons);
2020     evTag->SetNumOfElectrons(nElectrons);
2021     evTag->SetNumOfPhotons(nGamas);
2022     evTag->SetNumOfPi0s(nPi0s);
2023     evTag->SetNumOfNeutrons(nNeutrons);
2024     evTag->SetNumOfKaon0s(nK0s);
2025
2026     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
2027     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
2028     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
2029     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
2030     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
2031     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
2032     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
2033     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
2034     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
2035
2036     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
2037     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2038
2039     evTag->SetTotalMomentum(totalP);
2040     evTag->SetMeanPt(meanPt);
2041     evTag->SetMaxPt(maxPt);
2042
2043     tag->SetLHCTag(lhcLuminosity,lhcState);
2044     tag->SetDetectorTag(detectorMask);
2045
2046     tag->SetRunId(iInitRunNumber);
2047     tag->AddEventTag(*evTag);
2048   }
2049   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2050   else lastEvent = fLastEvent;
2051
2052   ttag.Fill();
2053   tag->Clear();
2054
2055   char fileName[256];
2056   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2057           tag->GetRunId(),fFirstEvent,lastEvent );
2058   AliInfo(Form("writing tags to file %s", fileName));
2059   AliDebug(1, Form("writing tags to file %s", fileName));
2060
2061   TFile* ftag = TFile::Open(fileName, "recreate");
2062   ftag->cd();
2063   ttag.Write();
2064   ftag->Close();
2065   file->cd();
2066   delete tag;
2067   delete evTag;
2068 }
2069
2070 //_____________________________________________________________________________
2071 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
2072 {
2073   // write all files from the given esd file to an aod file
2074
2075   // create an AliAOD object
2076   AliAODEvent *aod = new AliAODEvent();
2077   aod->CreateStdContent();
2078
2079   // go to the file
2080   aodFile->cd();
2081
2082   // create the tree
2083   TTree *aodTree = new TTree("AOD", "AliAOD tree");
2084   aodTree->Branch(aod->GetList());
2085
2086   // connect to ESD
2087   TTree *t = (TTree*) esdFile->Get("esdTree");
2088   TBranch *b = t->GetBranch("ESD");
2089   AliESD *esd = 0;
2090   b->SetAddress(&esd);
2091
2092   Int_t nEvents = b->GetEntries();
2093
2094   // set arrays and pointers
2095   Float_t posF[3];
2096   Double_t pos[3];
2097   Double_t p[3];
2098   Double_t covVtx[6];
2099   Double_t covTr[21];
2100   Double_t pid[10];
2101
2102   // loop over events and fill them
2103   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2104     b->GetEntry(iEvent);
2105
2106     // Multiplicity information needed by the header (to be revised!)
2107     Int_t nTracks   = esd->GetNumberOfTracks();
2108     Int_t nPosTracks = 0;
2109     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2110       if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2111
2112     // create the header
2113     aod->AddHeader(new AliAODHeader(esd->GetRunNumber(),
2114                                     esd->GetBunchCrossNumber(),
2115                                     esd->GetOrbitNumber(),
2116                                     esd->GetPeriodNumber(),
2117                                     nTracks,
2118                                     nPosTracks,
2119                                     nTracks-nPosTracks,
2120                                     esd->GetMagneticField(),
2121                                     -999., // fill muon magnetic field
2122                                     -999., // centrality; to be filled, still
2123                                     esd->GetZDCN1Energy(),
2124                                     esd->GetZDCP1Energy(),
2125                                     esd->GetZDCN2Energy(),
2126                                     esd->GetZDCP2Energy(),
2127                                     esd->GetZDCEMEnergy(),
2128                                     esd->GetTriggerMask(),
2129                                     esd->GetTriggerCluster(),
2130                                     esd->GetEventType()));
2131
2132     Int_t nV0s      = esd->GetNumberOfV0s();
2133     Int_t nCascades = esd->GetNumberOfCascades();
2134     Int_t nKinks    = esd->GetNumberOfKinks();
2135     Int_t nVertices = nV0s + nCascades + nKinks;
2136
2137     aod->ResetStd(nTracks, nVertices);
2138     AliAODTrack *aodTrack;
2139
2140
2141     // Array to take into account the tracks already added to the AOD
2142     Bool_t * usedTrack = NULL;
2143     if (nTracks>0) {
2144       usedTrack = new Bool_t[nTracks];
2145       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2146     }
2147     // Array to take into account the V0s already added to the AOD
2148     Bool_t * usedV0 = NULL;
2149     if (nV0s>0) {
2150       usedV0 = new Bool_t[nV0s];
2151       for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2152     }
2153     // Array to take into account the kinks already added to the AOD
2154     Bool_t * usedKink = NULL;
2155     if (nKinks>0) {
2156       usedKink = new Bool_t[nKinks];
2157       for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2158     }
2159
2160     // Access to the AOD container of vertices
2161     TClonesArray &vertices = *(aod->GetVertices());
2162     Int_t jVertices=0;
2163
2164     // Access to the AOD container of tracks
2165     TClonesArray &tracks = *(aod->GetTracks());
2166     Int_t jTracks=0;
2167
2168     // Add primary vertex. The primary tracks will be defined
2169     // after the loops on the composite objects (V0, cascades, kinks)
2170     const AliESDVertex *vtx = esd->GetPrimaryVertex();
2171
2172     vtx->GetXYZ(pos); // position
2173     vtx->GetCovMatrix(covVtx); //covariance matrix
2174
2175     AliAODVertex * primary = new(vertices[jVertices++])
2176       AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2177
2178     // Create vertices starting from the most complex objects
2179
2180     // Cascades
2181     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2182       AliESDcascade *cascade = esd->GetCascade(nCascade);
2183
2184       cascade->GetXYZ(pos[0], pos[1], pos[2]);
2185       cascade->GetPosCovXi(covVtx);
2186
2187       // Add the cascade vertex
2188       AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2189                                                                         covVtx,
2190                                                                         cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2191                                                                         primary,
2192                                                                         AliAODVertex::kCascade);
2193
2194       primary->AddDaughter(vcascade);
2195
2196       // Add the V0 from the cascade. The ESD class have to be optimized...
2197       // Now we have to search for the corresponding Vo in the list of V0s
2198       // using the indeces of the positive and negative tracks
2199
2200       Int_t posFromV0 = cascade->GetPindex();
2201       Int_t negFromV0 = cascade->GetNindex();
2202
2203
2204       AliESDv0 * v0 = 0x0;
2205       Int_t indV0 = -1;
2206
2207       for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2208
2209         v0 = esd->GetV0(iV0);
2210         Int_t posV0 = v0->GetPindex();
2211         Int_t negV0 = v0->GetNindex();
2212
2213         if (posV0==posFromV0 && negV0==negFromV0) {
2214           indV0 = iV0;
2215           break;
2216         }
2217       }
2218
2219       AliAODVertex * vV0FromCascade = 0x0;
2220
2221       if (indV0>-1 && !usedV0[indV0] ) {
2222
2223         // the V0 exists in the array of V0s and is not used
2224
2225         usedV0[indV0] = kTRUE;
2226
2227         v0->GetXYZ(pos[0], pos[1], pos[2]);
2228         v0->GetPosCov(covVtx);
2229
2230         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2231                                                                  covVtx,
2232                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2233                                                                  vcascade,
2234                                                                  AliAODVertex::kV0);
2235       } else {
2236
2237         // the V0 doesn't exist in the array of V0s or was used
2238         cerr << "Error: event " << iEvent << " cascade " << nCascade
2239              << " The V0 " << indV0
2240              << " doesn't exist in the array of V0s or was used!" << endl;
2241
2242         cascade->GetXYZ(pos[0], pos[1], pos[2]);
2243         cascade->GetPosCov(covVtx);
2244
2245         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2246                                                                  covVtx,
2247                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2248                                                                  vcascade,
2249                                                                  AliAODVertex::kV0);
2250         vcascade->AddDaughter(vV0FromCascade);
2251       }
2252
2253       // Add the positive tracks from the V0
2254
2255       if (! usedTrack[posFromV0]) {
2256
2257         usedTrack[posFromV0] = kTRUE;
2258
2259         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2260         esdTrack->GetPxPyPz(p);
2261         esdTrack->GetXYZ(pos);
2262         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2263         esdTrack->GetESDpid(pid);
2264
2265         vV0FromCascade->AddDaughter(aodTrack =
2266                                     new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2267                                            esdTrack->GetLabel(),
2268                                            p,
2269                                            kTRUE,
2270                                            pos,
2271                                            kFALSE,
2272                                            covTr,
2273                                            (Short_t)esdTrack->GetSign(),
2274                                            esdTrack->GetITSClusterMap(),
2275                                            pid,
2276                                            vV0FromCascade,
2277                                            kTRUE,  // check if this is right
2278                                            kFALSE, // check if this is right
2279                                            AliAODTrack::kSecondary)
2280                 );
2281         aodTrack->ConvertAliPIDtoAODPID();
2282       }
2283       else {
2284         cerr << "Error: event " << iEvent << " cascade " << nCascade
2285              << " track " << posFromV0 << " has already been used!" << endl;
2286       }
2287
2288       // Add the negative tracks from the V0
2289
2290       if (!usedTrack[negFromV0]) {
2291
2292         usedTrack[negFromV0] = kTRUE;
2293
2294         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2295         esdTrack->GetPxPyPz(p);
2296         esdTrack->GetXYZ(pos);
2297         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2298         esdTrack->GetESDpid(pid);
2299
2300         vV0FromCascade->AddDaughter(aodTrack =
2301                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2302                                            esdTrack->GetLabel(),
2303                                            p,
2304                                            kTRUE,
2305                                            pos,
2306                                            kFALSE,
2307                                            covTr,
2308                                            (Short_t)esdTrack->GetSign(),
2309                                            esdTrack->GetITSClusterMap(),
2310                                            pid,
2311                                            vV0FromCascade,
2312                                            kTRUE,  // check if this is right
2313                                            kFALSE, // check if this is right
2314                                            AliAODTrack::kSecondary)
2315                 );
2316         aodTrack->ConvertAliPIDtoAODPID();
2317       }
2318       else {
2319         cerr << "Error: event " << iEvent << " cascade " << nCascade
2320              << " track " << negFromV0 << " has already been used!" << endl;
2321       }
2322
2323       // Add the bachelor track from the cascade
2324
2325       Int_t bachelor = cascade->GetBindex();
2326
2327       if(!usedTrack[bachelor]) {
2328
2329         usedTrack[bachelor] = kTRUE;
2330
2331         AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2332         esdTrack->GetPxPyPz(p);
2333         esdTrack->GetXYZ(pos);
2334         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2335         esdTrack->GetESDpid(pid);
2336
2337         vcascade->AddDaughter(aodTrack =
2338                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2339                                            esdTrack->GetLabel(),
2340                                            p,
2341                                            kTRUE,
2342                                            pos,
2343                                            kFALSE,
2344                                            covTr,
2345                                            (Short_t)esdTrack->GetSign(),
2346                                            esdTrack->GetITSClusterMap(),
2347                                            pid,
2348                                            vcascade,
2349                                            kTRUE,  // check if this is right
2350                                            kFALSE, // check if this is right
2351                                            AliAODTrack::kSecondary)
2352                 );
2353         aodTrack->ConvertAliPIDtoAODPID();
2354      }
2355       else {
2356         cerr << "Error: event " << iEvent << " cascade " << nCascade
2357              << " track " << bachelor << " has already been used!" << endl;
2358       }
2359
2360       // Add the primary track of the cascade (if any)
2361
2362     } // end of the loop on cascades
2363
2364     // V0s
2365
2366     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2367
2368       if (usedV0[nV0]) continue; // skip if aready added to the AOD
2369
2370       AliESDv0 *v0 = esd->GetV0(nV0);
2371
2372       v0->GetXYZ(pos[0], pos[1], pos[2]);
2373       v0->GetPosCov(covVtx);
2374
2375       AliAODVertex * vV0 =
2376         new(vertices[jVertices++]) AliAODVertex(pos,
2377                                                 covVtx,
2378                                                 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2379                                                 primary,
2380                                                 AliAODVertex::kV0);
2381       primary->AddDaughter(vV0);
2382
2383       Int_t posFromV0 = v0->GetPindex();
2384       Int_t negFromV0 = v0->GetNindex();
2385
2386       // Add the positive tracks from the V0
2387
2388       if (!usedTrack[posFromV0]) {
2389
2390         usedTrack[posFromV0] = kTRUE;
2391
2392         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2393         esdTrack->GetPxPyPz(p);
2394         esdTrack->GetXYZ(pos);
2395         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2396         esdTrack->GetESDpid(pid);
2397
2398         vV0->AddDaughter(aodTrack =
2399                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2400                                            esdTrack->GetLabel(),
2401                                            p,
2402                                            kTRUE,
2403                                            pos,
2404                                            kFALSE,
2405                                            covTr,
2406                                            (Short_t)esdTrack->GetSign(),
2407                                            esdTrack->GetITSClusterMap(),
2408                                            pid,
2409                                            vV0,
2410                                            kTRUE,  // check if this is right
2411                                            kFALSE, // check if this is right
2412                                            AliAODTrack::kSecondary)
2413                 );
2414         aodTrack->ConvertAliPIDtoAODPID();
2415       }
2416       else {
2417         cerr << "Error: event " << iEvent << " V0 " << nV0
2418              << " track " << posFromV0 << " has already been used!" << endl;
2419       }
2420
2421       // Add the negative tracks from the V0
2422
2423       if (!usedTrack[negFromV0]) {
2424
2425         usedTrack[negFromV0] = kTRUE;
2426
2427         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2428         esdTrack->GetPxPyPz(p);
2429         esdTrack->GetXYZ(pos);
2430         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2431         esdTrack->GetESDpid(pid);
2432
2433         vV0->AddDaughter(aodTrack =
2434                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2435                                            esdTrack->GetLabel(),
2436                                            p,
2437                                            kTRUE,
2438                                            pos,
2439                                            kFALSE,
2440                                            covTr,
2441                                            (Short_t)esdTrack->GetSign(),
2442                                            esdTrack->GetITSClusterMap(),
2443                                            pid,
2444                                            vV0,
2445                                            kTRUE,  // check if this is right
2446                                            kFALSE, // check if this is right
2447                                            AliAODTrack::kSecondary)
2448                 );
2449         aodTrack->ConvertAliPIDtoAODPID();
2450       }
2451       else {
2452         cerr << "Error: event " << iEvent << " V0 " << nV0
2453              << " track " << negFromV0 << " has already been used!" << endl;
2454       }
2455
2456     } // end of the loop on V0s
2457
2458     // Kinks: it is a big mess the access to the information in the kinks
2459     // The loop is on the tracks in order to find the mother and daugther of each kink
2460
2461
2462     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2463
2464
2465       AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2466
2467       Int_t ikink = esdTrack->GetKinkIndex(0);
2468
2469       if (ikink) {
2470         // Negative kink index: mother, positive: daughter
2471
2472         // Search for the second track of the kink
2473
2474         for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2475
2476           AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2477
2478           Int_t jkink = esdTrack1->GetKinkIndex(0);
2479
2480           if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2481
2482             // The two tracks are from the same kink
2483
2484             if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2485
2486             Int_t imother = -1;
2487             Int_t idaughter = -1;
2488
2489             if (ikink<0 && jkink>0) {
2490
2491               imother = iTrack;
2492               idaughter = jTrack;
2493             }
2494             else if (ikink>0 && jkink<0) {
2495
2496               imother = jTrack;
2497               idaughter = iTrack;
2498             }
2499             else {
2500               cerr << "Error: Wrong combination of kink indexes: "
2501               << ikink << " " << jkink << endl;
2502               continue;
2503             }
2504
2505             // Add the mother track
2506
2507             AliAODTrack * mother = NULL;
2508
2509             if (!usedTrack[imother]) {
2510
2511               usedTrack[imother] = kTRUE;
2512
2513               AliESDtrack *esdTrack = esd->GetTrack(imother);
2514               esdTrack->GetPxPyPz(p);
2515               esdTrack->GetXYZ(pos);
2516               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2517               esdTrack->GetESDpid(pid);
2518
2519               mother =
2520                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2521                                            esdTrack->GetLabel(),
2522                                            p,
2523                                            kTRUE,
2524                                            pos,
2525                                            kFALSE,
2526                                            covTr,
2527                                            (Short_t)esdTrack->GetSign(),
2528                                            esdTrack->GetITSClusterMap(),
2529                                            pid,
2530                                            primary,
2531                                            kTRUE, // check if this is right
2532                                            kTRUE, // check if this is right
2533                                            AliAODTrack::kPrimary);
2534               primary->AddDaughter(mother);
2535               mother->ConvertAliPIDtoAODPID();
2536             }
2537             else {
2538               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2539               << " track " << imother << " has already been used!" << endl;
2540             }
2541
2542             // Add the kink vertex
2543             AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2544
2545             AliAODVertex * vkink =
2546             new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2547                                                     NULL,
2548                                                     0.,
2549                                                     mother,
2550                                                     AliAODVertex::kKink);
2551             // Add the daughter track
2552
2553             AliAODTrack * daughter = NULL;
2554
2555             if (!usedTrack[idaughter]) {
2556
2557               usedTrack[idaughter] = kTRUE;
2558
2559               AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2560               esdTrack->GetPxPyPz(p);
2561               esdTrack->GetXYZ(pos);
2562               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2563               esdTrack->GetESDpid(pid);
2564
2565               daughter =
2566                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2567                                            esdTrack->GetLabel(),
2568                                            p,
2569                                            kTRUE,
2570                                            pos,
2571                                            kFALSE,
2572                                            covTr,
2573                                            (Short_t)esdTrack->GetSign(),
2574                                            esdTrack->GetITSClusterMap(),
2575                                            pid,
2576                                            vkink,
2577                                            kTRUE, // check if this is right
2578                                            kTRUE, // check if this is right
2579                                            AliAODTrack::kPrimary);
2580               vkink->AddDaughter(daughter);
2581               daughter->ConvertAliPIDtoAODPID();
2582             }
2583             else {
2584               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2585               << " track " << idaughter << " has already been used!" << endl;
2586             }
2587
2588
2589           }
2590         }
2591
2592       }
2593
2594     }
2595
2596
2597     // Tracks (primary and orphan)
2598
2599     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2600
2601
2602       if (usedTrack[nTrack]) continue;
2603
2604       AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2605       esdTrack->GetPxPyPz(p);
2606       esdTrack->GetXYZ(pos);
2607       esdTrack->GetCovarianceXYZPxPyPz(covTr);
2608       esdTrack->GetESDpid(pid);
2609
2610       Float_t impactXY, impactZ;
2611
2612       esdTrack->GetImpactParameters(impactXY,impactZ);
2613
2614       if (impactXY<3) {
2615         // track inside the beam pipe
2616
2617         primary->AddDaughter(aodTrack =
2618             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2619                                          esdTrack->GetLabel(),
2620                                          p,
2621                                          kTRUE,
2622                                          pos,
2623                                          kFALSE,
2624                                          covTr,
2625                                          (Short_t)esdTrack->GetSign(),
2626                                          esdTrack->GetITSClusterMap(),
2627                                          pid,
2628                                          primary,
2629                                          kTRUE, // check if this is right
2630                                          kTRUE, // check if this is right
2631                                          AliAODTrack::kPrimary)
2632             );
2633         aodTrack->ConvertAliPIDtoAODPID();
2634       }
2635       else {
2636         // outside the beam pipe: orphan track
2637             aodTrack =
2638             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2639                                          esdTrack->GetLabel(),
2640                                          p,
2641                                          kTRUE,
2642                                          pos,
2643                                          kFALSE,
2644                                          covTr,
2645                                          (Short_t)esdTrack->GetSign(),
2646                                          esdTrack->GetITSClusterMap(),
2647                                          pid,
2648                                          NULL,
2649                                          kFALSE, // check if this is right
2650                                          kFALSE, // check if this is right
2651                                          AliAODTrack::kOrphan);
2652             aodTrack->ConvertAliPIDtoAODPID();
2653       }
2654     } // end of loop on tracks
2655
2656     // muon tracks
2657     Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2658     for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2659
2660       AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2661       p[0] = esdMuTrack->Px();
2662       p[1] = esdMuTrack->Py();
2663       p[2] = esdMuTrack->Pz();
2664       pos[0] = primary->GetX();
2665       pos[1] = primary->GetY();
2666       pos[2] = primary->GetZ();
2667
2668       // has to be changed once the muon pid is provided by the ESD
2669       for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2670
2671       primary->AddDaughter(
2672           new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2673                                              0, // no label provided
2674                                              p,
2675                                              kTRUE,
2676                                              pos,
2677                                              kFALSE,
2678                                              NULL, // no covariance matrix provided
2679                                              (Short_t)-99, // no charge provided
2680                                              0, // no ITSClusterMap
2681                                              pid,
2682                                              primary,
2683                                              kTRUE,  // check if this is right
2684                                              kTRUE,  // not used for vertex fit
2685                                              AliAODTrack::kPrimary)
2686           );
2687     }
2688
2689     // Access to the AOD container of clusters
2690     TClonesArray &clusters = *(aod->GetClusters());
2691     Int_t jClusters=0;
2692
2693     // Calo Clusters
2694     Int_t nClusters    = esd->GetNumberOfCaloClusters();
2695
2696     for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2697
2698       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2699
2700       Int_t id = cluster->GetID();
2701       Int_t label = -1;
2702       Float_t energy = cluster->GetClusterEnergy();
2703       cluster->GetGlobalPosition(posF);
2704       AliAODVertex *prodVertex = primary;
2705       AliAODTrack *primTrack = NULL;
2706       Char_t ttype=AliAODCluster::kUndef;
2707
2708       if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2709       else if (cluster->IsEMCAL()) {
2710
2711         if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2712           ttype = AliAODCluster::kEMCALPseudoCluster;
2713         else
2714           ttype = AliAODCluster::kEMCALClusterv1;
2715
2716       }
2717
2718       new(clusters[jClusters++]) AliAODCluster(id,
2719                                                label,
2720                                                energy,
2721                                                pos,
2722                                                NULL, // no covariance matrix provided
2723                                                NULL, // no pid for clusters provided
2724                                                prodVertex,
2725                                                primTrack,
2726                                                ttype);
2727
2728     } // end of loop on calo clusters
2729
2730     delete [] usedTrack;
2731     delete [] usedV0;
2732     delete [] usedKink;
2733
2734     // fill the tree for this event
2735     aodTree->Fill();
2736   } // end of event loop
2737
2738   aodTree->GetUserInfo()->Add(aod);
2739
2740   // close ESD file
2741   esdFile->Close();
2742
2743   // write the tree to the specified file
2744   aodFile = aodTree->GetCurrentFile();
2745   aodFile->cd();
2746   aodTree->Write();
2747
2748   return;
2749 }
2750
2751
2752 void AliReconstruction::WriteAlignmentData(AliESD* esd)
2753 {
2754   // Write space-points which are then used in the alignment procedures
2755   // For the moment only ITS, TRD and TPC
2756
2757   // Load TOF clusters
2758   if (fTracker[3]){
2759     fLoader[3]->LoadRecPoints("read");
2760     TTree* tree = fLoader[3]->TreeR();
2761     if (!tree) {
2762       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2763       return;
2764     }
2765     fTracker[3]->LoadClusters(tree);
2766   }
2767   Int_t ntracks = esd->GetNumberOfTracks();
2768   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2769     {
2770       AliESDtrack *track = esd->GetTrack(itrack);
2771       Int_t nsp = 0;
2772       Int_t idx[200];
2773       for (Int_t iDet = 3; iDet >= 0; iDet--)
2774         nsp += track->GetNcls(iDet);
2775       if (nsp) {
2776         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2777         track->SetTrackPointArray(sp);
2778         Int_t isptrack = 0;
2779         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2780           AliTracker *tracker = fTracker[iDet];
2781           if (!tracker) continue;
2782           Int_t nspdet = track->GetNcls(iDet);
2783           if (nspdet <= 0) continue;
2784           track->GetClusters(iDet,idx);
2785           AliTrackPoint p;
2786           Int_t isp = 0;
2787           Int_t isp2 = 0;
2788           while (isp < nspdet) {
2789             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2790             const Int_t kNTPCmax = 159;
2791             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2792             if (!isvalid) continue;
2793             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2794           }
2795         }
2796       }
2797     }
2798   if (fTracker[3]){
2799     fTracker[3]->UnloadClusters();
2800     fLoader[3]->UnloadRecPoints();
2801   }
2802 }
2803
2804 //_____________________________________________________________________________
2805 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2806 {
2807   // The method reads the raw-data error log
2808   // accumulated within the rawReader.
2809   // It extracts the raw-data errors related to
2810   // the current event and stores them into
2811   // a TClonesArray inside the esd object.
2812
2813   if (!fRawReader) return;
2814
2815   for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2816
2817     AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2818     if (!log) continue;
2819     if (iEvent != log->GetEventNumber()) continue;
2820
2821     esd->AddRawDataErrorLog(log);
2822   }
2823
2824 }