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   if (fVertexer) {
1045   AliMultiplicity *mult= fVertexer->GetMultiplicity();
1046   if(mult)esd->SetMultiplicity(mult);
1047   }
1048
1049   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1050     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1051   }
1052   delete vertex;
1053
1054   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1055                stopwatch.RealTime(),stopwatch.CpuTime()));
1056
1057   return kTRUE;
1058 }
1059
1060 //_____________________________________________________________________________
1061 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1062 {
1063 // run the HLT barrel tracking
1064
1065   TStopwatch stopwatch;
1066   stopwatch.Start();
1067
1068   if (!fRunLoader) {
1069     AliError("Missing runLoader!");
1070     return kFALSE;
1071   }
1072
1073   AliInfo("running HLT tracking");
1074
1075   // Get a pointer to the HLT reconstructor
1076   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1077   if (!reconstructor) return kFALSE;
1078
1079   // TPC + ITS
1080   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1081     TString detName = fgkDetectorName[iDet];
1082     AliDebug(1, Form("%s HLT tracking", detName.Data()));
1083     reconstructor->SetOption(detName.Data());
1084     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1085     if (!tracker) {
1086       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1087       if (fStopOnError) return kFALSE;
1088       continue;
1089     }
1090     Double_t vtxPos[3];
1091     Double_t vtxErr[3]={0.005,0.005,0.010};
1092     const AliESDVertex *vertex = esd->GetVertex();
1093     vertex->GetXYZ(vtxPos);
1094     tracker->SetVertex(vtxPos,vtxErr);
1095     if(iDet != 1) {
1096       fLoader[iDet]->LoadRecPoints("read");
1097       TTree* tree = fLoader[iDet]->TreeR();
1098       if (!tree) {
1099         AliError(Form("Can't get the %s cluster tree", detName.Data()));
1100         return kFALSE;
1101       }
1102       tracker->LoadClusters(tree);
1103     }
1104     if (tracker->Clusters2Tracks(esd) != 0) {
1105       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1106       return kFALSE;
1107     }
1108     if(iDet != 1) {
1109       tracker->UnloadClusters();
1110     }
1111     delete tracker;
1112   }
1113
1114   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1115                stopwatch.RealTime(),stopwatch.CpuTime()));
1116
1117   return kTRUE;
1118 }
1119
1120 //_____________________________________________________________________________
1121 Bool_t AliReconstruction::RunMuonTracking()
1122 {
1123 // run the muon spectrometer tracking
1124
1125   TStopwatch stopwatch;
1126   stopwatch.Start();
1127
1128   if (!fRunLoader) {
1129     AliError("Missing runLoader!");
1130     return kFALSE;
1131   }
1132   Int_t iDet = 7; // for MUON
1133
1134   AliInfo("is running...");
1135
1136   // Get a pointer to the MUON reconstructor
1137   AliReconstructor *reconstructor = GetReconstructor(iDet);
1138   if (!reconstructor) return kFALSE;
1139
1140
1141   TString detName = fgkDetectorName[iDet];
1142   AliDebug(1, Form("%s tracking", detName.Data()));
1143   AliTracker *tracker =  reconstructor->CreateTracker(fRunLoader);
1144   if (!tracker) {
1145     AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1146     return kFALSE;
1147   }
1148
1149   // create Tracks
1150   fLoader[iDet]->LoadTracks("update");
1151   fLoader[iDet]->CleanTracks();
1152   fLoader[iDet]->MakeTracksContainer();
1153
1154   // read RecPoints
1155   fLoader[iDet]->LoadRecPoints("read");
1156
1157   if (!tracker->Clusters2Tracks(0x0)) {
1158     AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1159     return kFALSE;
1160   }
1161   fLoader[iDet]->UnloadRecPoints();
1162
1163   fLoader[iDet]->WriteTracks("OVERWRITE");
1164   fLoader[iDet]->UnloadTracks();
1165
1166   delete tracker;
1167
1168
1169   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1170                stopwatch.RealTime(),stopwatch.CpuTime()));
1171
1172   return kTRUE;
1173 }
1174
1175
1176 //_____________________________________________________________________________
1177 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1178 {
1179 // run the barrel tracking
1180
1181   TStopwatch stopwatch;
1182   stopwatch.Start();
1183
1184   AliInfo("running tracking");
1185
1186   //Fill the ESD with the T0 info (will be used by the TOF)
1187   if (fReconstructor[11])
1188       GetReconstructor(11)->FillESD(fRunLoader, esd);
1189
1190   // pass 1: TPC + ITS inwards
1191   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1192     if (!fTracker[iDet]) continue;
1193     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1194
1195     // load clusters
1196     fLoader[iDet]->LoadRecPoints("read");
1197     TTree* tree = fLoader[iDet]->TreeR();
1198     if (!tree) {
1199       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1200       return kFALSE;
1201     }
1202     fTracker[iDet]->LoadClusters(tree);
1203
1204     // run tracking
1205     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1206       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1207       return kFALSE;
1208     }
1209     if (fCheckPointLevel > 1) {
1210       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1211     }
1212     // preliminary PID in TPC needed by the ITS tracker
1213     if (iDet == 1) {
1214       GetReconstructor(1)->FillESD(fRunLoader, esd);
1215       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1216       AliESDpid::MakePID(esd);
1217     }
1218   }
1219
1220   // pass 2: ALL backwards
1221   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1222     if (!fTracker[iDet]) continue;
1223     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1224
1225     // load clusters
1226     if (iDet > 1) {     // all except ITS, TPC
1227       TTree* tree = NULL;
1228       fLoader[iDet]->LoadRecPoints("read");
1229       tree = fLoader[iDet]->TreeR();
1230       if (!tree) {
1231         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1232         return kFALSE;
1233       }
1234       fTracker[iDet]->LoadClusters(tree);
1235     }
1236
1237     // run tracking
1238     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1239       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1240       return kFALSE;
1241     }
1242     if (fCheckPointLevel > 1) {
1243       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1244     }
1245
1246     // unload clusters
1247     if (iDet > 2) {     // all except ITS, TPC, TRD
1248       fTracker[iDet]->UnloadClusters();
1249       fLoader[iDet]->UnloadRecPoints();
1250     }
1251     // updated PID in TPC needed by the ITS tracker -MI
1252     if (iDet == 1) {
1253       GetReconstructor(1)->FillESD(fRunLoader, esd);
1254       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1255       AliESDpid::MakePID(esd);
1256     }
1257   }
1258
1259   // write space-points to the ESD in case alignment data output
1260   // is switched on
1261   if (fWriteAlignmentData)
1262     WriteAlignmentData(esd);
1263
1264   // pass 3: TRD + TPC + ITS refit inwards
1265   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1266     if (!fTracker[iDet]) continue;
1267     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1268
1269     // run tracking
1270     if (fTracker[iDet]->RefitInward(esd) != 0) {
1271       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1272       return kFALSE;
1273     }
1274     if (fCheckPointLevel > 1) {
1275       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1276     }
1277
1278     // unload clusters
1279     fTracker[iDet]->UnloadClusters();
1280     fLoader[iDet]->UnloadRecPoints();
1281   }
1282   //
1283   // Propagate track to the vertex - if not done by ITS
1284   //
1285   Int_t ntracks = esd->GetNumberOfTracks();
1286   for (Int_t itrack=0; itrack<ntracks; itrack++){
1287     const Double_t kRadius  = 3;   // beam pipe radius
1288     const Double_t kMaxStep = 5;   // max step
1289     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1290     Double_t       fieldZ   = AliTracker::GetBz();  //
1291     AliESDtrack * track = esd->GetTrack(itrack);
1292     if (!track) continue;
1293     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1294     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1295     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1296   }
1297
1298   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1299                stopwatch.RealTime(),stopwatch.CpuTime()));
1300
1301   return kTRUE;
1302 }
1303
1304 //_____________________________________________________________________________
1305 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1306 {
1307 // fill the event summary data
1308
1309   TStopwatch stopwatch;
1310   stopwatch.Start();
1311   AliInfo("filling ESD");
1312
1313   TString detStr = detectors;
1314   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1315     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1316     AliReconstructor* reconstructor = GetReconstructor(iDet);
1317     if (!reconstructor) continue;
1318
1319     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1320       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1321       TTree* clustersTree = NULL;
1322       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1323         fLoader[iDet]->LoadRecPoints("read");
1324         clustersTree = fLoader[iDet]->TreeR();
1325         if (!clustersTree) {
1326           AliError(Form("Can't get the %s clusters tree",
1327                         fgkDetectorName[iDet]));
1328           if (fStopOnError) return kFALSE;
1329         }
1330       }
1331       if (fRawReader && !reconstructor->HasDigitConversion()) {
1332         reconstructor->FillESD(fRawReader, clustersTree, esd);
1333       } else {
1334         TTree* digitsTree = NULL;
1335         if (fLoader[iDet]) {
1336           fLoader[iDet]->LoadDigits("read");
1337           digitsTree = fLoader[iDet]->TreeD();
1338           if (!digitsTree) {
1339             AliError(Form("Can't get the %s digits tree",
1340                           fgkDetectorName[iDet]));
1341             if (fStopOnError) return kFALSE;
1342           }
1343         }
1344         reconstructor->FillESD(digitsTree, clustersTree, esd);
1345         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1346       }
1347       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1348         fLoader[iDet]->UnloadRecPoints();
1349       }
1350
1351       if (fRawReader) {
1352         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1353       } else {
1354         reconstructor->FillESD(fRunLoader, esd);
1355       }
1356       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1357     }
1358   }
1359
1360   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1361     AliError(Form("the following detectors were not found: %s",
1362                   detStr.Data()));
1363     if (fStopOnError) return kFALSE;
1364   }
1365
1366   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1367                stopwatch.RealTime(),stopwatch.CpuTime()));
1368
1369   return kTRUE;
1370 }
1371
1372 //_____________________________________________________________________________
1373 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1374 {
1375   // Reads the trigger decision which is
1376   // stored in Trigger.root file and fills
1377   // the corresponding esd entries
1378
1379   AliInfo("Filling trigger information into the ESD");
1380
1381   if (fRawReader) {
1382     AliCTPRawStream input(fRawReader);
1383     if (!input.Next()) {
1384       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1385       return kFALSE;
1386     }
1387     esd->SetTriggerMask(input.GetClassMask());
1388     esd->SetTriggerCluster(input.GetClusterMask());
1389   }
1390   else {
1391     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1392     if (runloader) {
1393       if (!runloader->LoadTrigger()) {
1394         AliCentralTrigger *aCTP = runloader->GetTrigger();
1395         esd->SetTriggerMask(aCTP->GetClassMask());
1396         esd->SetTriggerCluster(aCTP->GetClusterMask());
1397       }
1398       else {
1399         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1400         return kFALSE;
1401       }
1402     }
1403     else {
1404       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1405       return kFALSE;
1406     }
1407   }
1408
1409   return kTRUE;
1410 }
1411
1412
1413
1414
1415
1416 //_____________________________________________________________________________
1417 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1418 {
1419   //
1420   // Filling information from RawReader Header
1421   //
1422
1423   AliInfo("Filling information from RawReader Header");
1424   esd->SetBunchCrossNumber(0);
1425   esd->SetOrbitNumber(0);
1426   esd->SetPeriodNumber(0);
1427   esd->SetTimeStamp(0);
1428   esd->SetEventType(0);
1429   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1430   if (eventHeader){
1431
1432     const UInt_t *id = eventHeader->GetP("Id");
1433     esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1434     esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1435     esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1436
1437     esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1438     esd->SetEventType((eventHeader->Get("Type")));
1439   }
1440
1441   return kTRUE;
1442 }
1443
1444
1445 //_____________________________________________________________________________
1446 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1447 {
1448 // check whether detName is contained in detectors
1449 // if yes, it is removed from detectors
1450
1451   // check if all detectors are selected
1452   if ((detectors.CompareTo("ALL") == 0) ||
1453       detectors.BeginsWith("ALL ") ||
1454       detectors.EndsWith(" ALL") ||
1455       detectors.Contains(" ALL ")) {
1456     detectors = "ALL";
1457     return kTRUE;
1458   }
1459
1460   // search for the given detector
1461   Bool_t result = kFALSE;
1462   if ((detectors.CompareTo(detName) == 0) ||
1463       detectors.BeginsWith(detName+" ") ||
1464       detectors.EndsWith(" "+detName) ||
1465       detectors.Contains(" "+detName+" ")) {
1466     detectors.ReplaceAll(detName, "");
1467     result = kTRUE;
1468   }
1469
1470   // clean up the detectors string
1471   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1472   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1473   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1474
1475   return result;
1476 }
1477
1478 //_____________________________________________________________________________
1479 Bool_t AliReconstruction::InitRunLoader()
1480 {
1481 // get or create the run loader
1482
1483   if (gAlice) delete gAlice;
1484   gAlice = NULL;
1485
1486   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1487     // load all base libraries to get the loader classes
1488     TString libs = gSystem->GetLibraries();
1489     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1490       TString detName = fgkDetectorName[iDet];
1491       if (detName == "HLT") continue;
1492       if (libs.Contains("lib" + detName + "base.so")) continue;
1493       gSystem->Load("lib" + detName + "base.so");
1494     }
1495     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1496     if (!fRunLoader) {
1497       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1498       CleanUp();
1499       return kFALSE;
1500     }
1501     fRunLoader->CdGAFile();
1502     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1503       if (fRunLoader->LoadgAlice() == 0) {
1504         gAlice = fRunLoader->GetAliRun();
1505         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1506       }
1507     }
1508     if (!gAlice && !fRawReader) {
1509       AliError(Form("no gAlice object found in file %s",
1510                     fGAliceFileName.Data()));
1511       CleanUp();
1512       return kFALSE;
1513     }
1514
1515   } else {               // galice.root does not exist
1516     if (!fRawReader) {
1517       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1518       CleanUp();
1519       return kFALSE;
1520     }
1521     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1522                                     AliConfig::GetDefaultEventFolderName(),
1523                                     "recreate");
1524     if (!fRunLoader) {
1525       AliError(Form("could not create run loader in file %s",
1526                     fGAliceFileName.Data()));
1527       CleanUp();
1528       return kFALSE;
1529     }
1530     fRunLoader->MakeTree("E");
1531     Int_t iEvent = 0;
1532     while (fRawReader->NextEvent()) {
1533       fRunLoader->SetEventNumber(iEvent);
1534       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1535                                      iEvent, iEvent);
1536       fRunLoader->MakeTree("H");
1537       fRunLoader->TreeE()->Fill();
1538       iEvent++;
1539     }
1540     fRawReader->RewindEvents();
1541     if (fNumberOfEventsPerFile > 0)
1542       fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1543     else
1544       fRunLoader->SetNumberOfEventsPerFile(iEvent);
1545     fRunLoader->WriteHeader("OVERWRITE");
1546     fRunLoader->CdGAFile();
1547     fRunLoader->Write(0, TObject::kOverwrite);
1548 //    AliTracker::SetFieldMap(???);
1549   }
1550
1551   return kTRUE;
1552 }
1553
1554 //_____________________________________________________________________________
1555 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1556 {
1557 // get the reconstructor object and the loader for a detector
1558
1559   if (fReconstructor[iDet]) return fReconstructor[iDet];
1560
1561   // load the reconstructor object
1562   TPluginManager* pluginManager = gROOT->GetPluginManager();
1563   TString detName = fgkDetectorName[iDet];
1564   TString recName = "Ali" + detName + "Reconstructor";
1565   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1566
1567   if (detName == "HLT") {
1568     if (!gROOT->GetClass("AliLevel3")) {
1569       gSystem->Load("libAliHLTSrc.so");
1570       gSystem->Load("libAliHLTMisc.so");
1571       gSystem->Load("libAliHLTHough.so");
1572       gSystem->Load("libAliHLTComp.so");
1573     }
1574   }
1575
1576   AliReconstructor* reconstructor = NULL;
1577   // first check if a plugin is defined for the reconstructor
1578   TPluginHandler* pluginHandler =
1579     pluginManager->FindHandler("AliReconstructor", detName);
1580   // if not, add a plugin for it
1581   if (!pluginHandler) {
1582     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1583     TString libs = gSystem->GetLibraries();
1584     if (libs.Contains("lib" + detName + "base.so") ||
1585         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1586       pluginManager->AddHandler("AliReconstructor", detName,
1587                                 recName, detName + "rec", recName + "()");
1588     } else {
1589       pluginManager->AddHandler("AliReconstructor", detName,
1590                                 recName, detName, recName + "()");
1591     }
1592     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1593   }
1594   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1595     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1596   }
1597   if (reconstructor) {
1598     TObject* obj = fOptions.FindObject(detName.Data());
1599     if (obj) reconstructor->SetOption(obj->GetTitle());
1600     reconstructor->Init(fRunLoader);
1601     fReconstructor[iDet] = reconstructor;
1602   }
1603
1604   // get or create the loader
1605   if (detName != "HLT") {
1606     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1607     if (!fLoader[iDet]) {
1608       AliConfig::Instance()
1609         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1610                                 detName, detName);
1611       // first check if a plugin is defined for the loader
1612       pluginHandler =
1613         pluginManager->FindHandler("AliLoader", detName);
1614       // if not, add a plugin for it
1615       if (!pluginHandler) {
1616         TString loaderName = "Ali" + detName + "Loader";
1617         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1618         pluginManager->AddHandler("AliLoader", detName,
1619                                   loaderName, detName + "base",
1620                                   loaderName + "(const char*, TFolder*)");
1621         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1622       }
1623       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1624         fLoader[iDet] =
1625           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1626                                                  fRunLoader->GetEventFolder());
1627       }
1628       if (!fLoader[iDet]) {   // use default loader
1629         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1630       }
1631       if (!fLoader[iDet]) {
1632         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1633         if (fStopOnError) return NULL;
1634       } else {
1635         fRunLoader->AddLoader(fLoader[iDet]);
1636         fRunLoader->CdGAFile();
1637         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1638         fRunLoader->Write(0, TObject::kOverwrite);
1639       }
1640     }
1641   }
1642
1643   return reconstructor;
1644 }
1645
1646 //_____________________________________________________________________________
1647 Bool_t AliReconstruction::CreateVertexer()
1648 {
1649 // create the vertexer
1650
1651   fVertexer = NULL;
1652   AliReconstructor* itsReconstructor = GetReconstructor(0);
1653   if (itsReconstructor) {
1654     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1655   }
1656   if (!fVertexer) {
1657     AliWarning("couldn't create a vertexer for ITS");
1658     if (fStopOnError) return kFALSE;
1659   }
1660
1661   return kTRUE;
1662 }
1663
1664 //_____________________________________________________________________________
1665 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1666 {
1667 // create the trackers
1668
1669   TString detStr = detectors;
1670   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1671     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1672     AliReconstructor* reconstructor = GetReconstructor(iDet);
1673     if (!reconstructor) continue;
1674     TString detName = fgkDetectorName[iDet];
1675     if (detName == "HLT") {
1676       fRunHLTTracking = kTRUE;
1677       continue;
1678     }
1679     if (detName == "MUON") {
1680       fRunMuonTracking = kTRUE;
1681       continue;
1682     }
1683
1684
1685     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1686     if (!fTracker[iDet] && (iDet < 7)) {
1687       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1688       if (fStopOnError) return kFALSE;
1689     }
1690   }
1691
1692   return kTRUE;
1693 }
1694
1695 //_____________________________________________________________________________
1696 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1697 {
1698 // delete trackers and the run loader and close and delete the file
1699
1700   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1701     delete fReconstructor[iDet];
1702     fReconstructor[iDet] = NULL;
1703     fLoader[iDet] = NULL;
1704     delete fTracker[iDet];
1705     fTracker[iDet] = NULL;
1706   }
1707   delete fVertexer;
1708   fVertexer = NULL;
1709   delete fDiamondProfile;
1710   fDiamondProfile = NULL;
1711
1712   delete fRunLoader;
1713   fRunLoader = NULL;
1714   delete fRawReader;
1715   fRawReader = NULL;
1716
1717   if (file) {
1718     file->Close();
1719     delete file;
1720   }
1721
1722   if (fileOld) {
1723     fileOld->Close();
1724     delete fileOld;
1725     gSystem->Unlink("AliESDs.old.root");
1726   }
1727 }
1728
1729
1730 //_____________________________________________________________________________
1731 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1732 {
1733 // read the ESD event from a file
1734
1735   if (!esd) return kFALSE;
1736   char fileName[256];
1737   sprintf(fileName, "ESD_%d.%d_%s.root",
1738           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1739   if (gSystem->AccessPathName(fileName)) return kFALSE;
1740
1741   AliInfo(Form("reading ESD from file %s", fileName));
1742   AliDebug(1, Form("reading ESD from file %s", fileName));
1743   TFile* file = TFile::Open(fileName);
1744   if (!file || !file->IsOpen()) {
1745     AliError(Form("opening %s failed", fileName));
1746     delete file;
1747     return kFALSE;
1748   }
1749
1750   gROOT->cd();
1751   delete esd;
1752   esd = (AliESD*) file->Get("ESD");
1753   file->Close();
1754   delete file;
1755   return kTRUE;
1756 }
1757
1758 //_____________________________________________________________________________
1759 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1760 {
1761 // write the ESD event to a file
1762
1763   if (!esd) return;
1764   char fileName[256];
1765   sprintf(fileName, "ESD_%d.%d_%s.root",
1766           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1767
1768   AliDebug(1, Form("writing ESD to file %s", fileName));
1769   TFile* file = TFile::Open(fileName, "recreate");
1770   if (!file || !file->IsOpen()) {
1771     AliError(Form("opening %s failed", fileName));
1772   } else {
1773     esd->Write("ESD");
1774     file->Close();
1775   }
1776   delete file;
1777 }
1778
1779
1780
1781
1782 //_____________________________________________________________________________
1783 void AliReconstruction::CreateTag(TFile* file)
1784 {
1785   //GRP
1786   Float_t lhcLuminosity = 0.0;
1787   TString lhcState = "test";
1788   UInt_t detectorMask = 0;
1789
1790   /////////////
1791   //muon code//
1792   ////////////
1793   Double_t fMUONMASS = 0.105658369;
1794   //Variables
1795   Double_t fX,fY,fZ ;
1796   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1797   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1798   Int_t fCharge;
1799   TLorentzVector fEPvector;
1800
1801   Float_t fZVertexCut = 10.0;
1802   Float_t fRhoVertexCut = 2.0;
1803
1804   Float_t fLowPtCut = 1.0;
1805   Float_t fHighPtCut = 3.0;
1806   Float_t fVeryHighPtCut = 10.0;
1807   ////////////
1808
1809   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1810
1811   // Creates the tags for all the events in a given ESD file
1812   Int_t ntrack;
1813   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1814   Int_t nPos, nNeg, nNeutr;
1815   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1816   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1817   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1818   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1819   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1820   Int_t fVertexflag;
1821   Int_t iRunNumber = 0;
1822   TString fVertexName("default");
1823
1824   AliRunTag *tag = new AliRunTag();
1825   AliEventTag *evTag = new AliEventTag();
1826   TTree ttag("T","A Tree with event tags");
1827   TBranch * btag = ttag.Branch("AliTAG", &tag);
1828   btag->SetCompressionLevel(9);
1829
1830   AliInfo(Form("Creating the tags......."));
1831
1832   if (!file || !file->IsOpen()) {
1833     AliError(Form("opening failed"));
1834     delete file;
1835     return ;
1836   }
1837   Int_t lastEvent = 0;
1838   TTree *t = (TTree*) file->Get("esdTree");
1839   TBranch * b = t->GetBranch("ESD");
1840   AliESD *esd = 0;
1841   b->SetAddress(&esd);
1842
1843   b->GetEntry(fFirstEvent);
1844   Int_t iInitRunNumber = esd->GetRunNumber();
1845
1846   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1847   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1848   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1849     ntrack = 0;
1850     nPos = 0;
1851     nNeg = 0;
1852     nNeutr =0;
1853     nK0s = 0;
1854     nNeutrons = 0;
1855     nPi0s = 0;
1856     nGamas = 0;
1857     nProtons = 0;
1858     nKaons = 0;
1859     nPions = 0;
1860     nMuons = 0;
1861     nElectrons = 0;
1862     nCh1GeV = 0;
1863     nCh3GeV = 0;
1864     nCh10GeV = 0;
1865     nMu1GeV = 0;
1866     nMu3GeV = 0;
1867     nMu10GeV = 0;
1868     nEl1GeV = 0;
1869     nEl3GeV = 0;
1870     nEl10GeV = 0;
1871     maxPt = .0;
1872     meanPt = .0;
1873     totalP = .0;
1874     fVertexflag = 0;
1875
1876     b->GetEntry(iEventNumber);
1877     iRunNumber = esd->GetRunNumber();
1878     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1879     const AliESDVertex * vertexIn = esd->GetVertex();
1880     if (!vertexIn) AliError("ESD has not defined vertex.");
1881     if (vertexIn) fVertexName = vertexIn->GetName();
1882     if(fVertexName != "default") fVertexflag = 1;
1883     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1884       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1885       UInt_t status = esdTrack->GetStatus();
1886
1887       //select only tracks with ITS refit
1888       if ((status&AliESDtrack::kITSrefit)==0) continue;
1889       //select only tracks with TPC refit
1890       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1891
1892       //select only tracks with the "combined PID"
1893       if ((status&AliESDtrack::kESDpid)==0) continue;
1894       Double_t p[3];
1895       esdTrack->GetPxPyPz(p);
1896       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1897       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1898       totalP += momentum;
1899       meanPt += fPt;
1900       if(fPt > maxPt) maxPt = fPt;
1901
1902       if(esdTrack->GetSign() > 0) {
1903         nPos++;
1904         if(fPt > fLowPtCut) nCh1GeV++;
1905         if(fPt > fHighPtCut) nCh3GeV++;
1906         if(fPt > fVeryHighPtCut) nCh10GeV++;
1907       }
1908       if(esdTrack->GetSign() < 0) {
1909         nNeg++;
1910         if(fPt > fLowPtCut) nCh1GeV++;
1911         if(fPt > fHighPtCut) nCh3GeV++;
1912         if(fPt > fVeryHighPtCut) nCh10GeV++;
1913       }
1914       if(esdTrack->GetSign() == 0) nNeutr++;
1915
1916       //PID
1917       Double_t prob[5];
1918       esdTrack->GetESDpid(prob);
1919
1920       Double_t rcc = 0.0;
1921       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1922       if(rcc == 0.0) continue;
1923       //Bayes' formula
1924       Double_t w[5];
1925       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1926
1927       //protons
1928       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1929       //kaons
1930       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1931       //pions
1932       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1933       //electrons
1934       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1935         nElectrons++;
1936         if(fPt > fLowPtCut) nEl1GeV++;
1937         if(fPt > fHighPtCut) nEl3GeV++;
1938         if(fPt > fVeryHighPtCut) nEl10GeV++;
1939       }
1940       ntrack++;
1941     }//track loop
1942
1943     /////////////
1944     //muon code//
1945     ////////////
1946     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1947     // loop over all reconstructed tracks (also first track of combination)
1948     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1949       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1950       if (muonTrack == 0x0) continue;
1951
1952       // Coordinates at vertex
1953       fZ = muonTrack->GetZ();
1954       fY = muonTrack->GetBendingCoor();
1955       fX = muonTrack->GetNonBendingCoor();
1956
1957       fThetaX = muonTrack->GetThetaX();
1958       fThetaY = muonTrack->GetThetaY();
1959
1960       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1961       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1962       fPxRec = fPzRec * TMath::Tan(fThetaX);
1963       fPyRec = fPzRec * TMath::Tan(fThetaY);
1964       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1965
1966       //ChiSquare of the track if needed
1967       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1968       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1969       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1970
1971       // total number of muons inside a vertex cut
1972       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1973         nMuons++;
1974         if(fEPvector.Pt() > fLowPtCut) {
1975           nMu1GeV++;
1976           if(fEPvector.Pt() > fHighPtCut) {
1977             nMu3GeV++;
1978             if (fEPvector.Pt() > fVeryHighPtCut) {
1979               nMu10GeV++;
1980             }
1981           }
1982         }
1983       }
1984     }//muon track loop
1985
1986     // Fill the event tags
1987     if(ntrack != 0)
1988       meanPt = meanPt/ntrack;
1989
1990     evTag->SetEventId(iEventNumber+1);
1991     if (vertexIn) {
1992       evTag->SetVertexX(vertexIn->GetXv());
1993       evTag->SetVertexY(vertexIn->GetYv());
1994       evTag->SetVertexZ(vertexIn->GetZv());
1995       evTag->SetVertexZError(vertexIn->GetZRes());
1996     }
1997     evTag->SetVertexFlag(fVertexflag);
1998
1999     evTag->SetT0VertexZ(esd->GetT0zVertex());
2000
2001     evTag->SetTriggerMask(esd->GetTriggerMask());
2002     evTag->SetTriggerCluster(esd->GetTriggerCluster());
2003
2004     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
2005     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
2006     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
2007     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2008     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
2009     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
2010
2011
2012     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
2013     evTag->SetNumOfPosTracks(nPos);
2014     evTag->SetNumOfNegTracks(nNeg);
2015     evTag->SetNumOfNeutrTracks(nNeutr);
2016
2017     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
2018     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
2019     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
2020     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
2021
2022     evTag->SetNumOfProtons(nProtons);
2023     evTag->SetNumOfKaons(nKaons);
2024     evTag->SetNumOfPions(nPions);
2025     evTag->SetNumOfMuons(nMuons);
2026     evTag->SetNumOfElectrons(nElectrons);
2027     evTag->SetNumOfPhotons(nGamas);
2028     evTag->SetNumOfPi0s(nPi0s);
2029     evTag->SetNumOfNeutrons(nNeutrons);
2030     evTag->SetNumOfKaon0s(nK0s);
2031
2032     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
2033     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
2034     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
2035     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
2036     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
2037     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
2038     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
2039     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
2040     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
2041
2042     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
2043     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2044
2045     evTag->SetTotalMomentum(totalP);
2046     evTag->SetMeanPt(meanPt);
2047     evTag->SetMaxPt(maxPt);
2048
2049     tag->SetLHCTag(lhcLuminosity,lhcState);
2050     tag->SetDetectorTag(detectorMask);
2051
2052     tag->SetRunId(iInitRunNumber);
2053     tag->AddEventTag(*evTag);
2054   }
2055   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2056   else lastEvent = fLastEvent;
2057
2058   ttag.Fill();
2059   tag->Clear();
2060
2061   char fileName[256];
2062   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2063           tag->GetRunId(),fFirstEvent,lastEvent );
2064   AliInfo(Form("writing tags to file %s", fileName));
2065   AliDebug(1, Form("writing tags to file %s", fileName));
2066
2067   TFile* ftag = TFile::Open(fileName, "recreate");
2068   ftag->cd();
2069   ttag.Write();
2070   ftag->Close();
2071   file->cd();
2072   delete tag;
2073   delete evTag;
2074 }
2075
2076 //_____________________________________________________________________________
2077 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
2078 {
2079   // write all files from the given esd file to an aod file
2080
2081   // create an AliAOD object
2082   AliAODEvent *aod = new AliAODEvent();
2083   aod->CreateStdContent();
2084
2085   // go to the file
2086   aodFile->cd();
2087
2088   // create the tree
2089   TTree *aodTree = new TTree("AOD", "AliAOD tree");
2090   aodTree->Branch(aod->GetList());
2091
2092   // connect to ESD
2093   TTree *t = (TTree*) esdFile->Get("esdTree");
2094   TBranch *b = t->GetBranch("ESD");
2095   AliESD *esd = 0;
2096   b->SetAddress(&esd);
2097
2098   Int_t nEvents = b->GetEntries();
2099
2100   // set arrays and pointers
2101   Float_t posF[3];
2102   Double_t pos[3];
2103   Double_t p[3];
2104   Double_t covVtx[6];
2105   Double_t covTr[21];
2106   Double_t pid[10];
2107
2108   // loop over events and fill them
2109   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2110     b->GetEntry(iEvent);
2111
2112     // Multiplicity information needed by the header (to be revised!)
2113     Int_t nTracks   = esd->GetNumberOfTracks();
2114     Int_t nPosTracks = 0;
2115     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2116       if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2117
2118     // create the header
2119     aod->AddHeader(new AliAODHeader(esd->GetRunNumber(),
2120                                     esd->GetBunchCrossNumber(),
2121                                     esd->GetOrbitNumber(),
2122                                     esd->GetPeriodNumber(),
2123                                     nTracks,
2124                                     nPosTracks,
2125                                     nTracks-nPosTracks,
2126                                     esd->GetMagneticField(),
2127                                     -999., // fill muon magnetic field
2128                                     -999., // centrality; to be filled, still
2129                                     esd->GetZDCN1Energy(),
2130                                     esd->GetZDCP1Energy(),
2131                                     esd->GetZDCN2Energy(),
2132                                     esd->GetZDCP2Energy(),
2133                                     esd->GetZDCEMEnergy(),
2134                                     esd->GetTriggerMask(),
2135                                     esd->GetTriggerCluster(),
2136                                     esd->GetEventType()));
2137
2138     Int_t nV0s      = esd->GetNumberOfV0s();
2139     Int_t nCascades = esd->GetNumberOfCascades();
2140     Int_t nKinks    = esd->GetNumberOfKinks();
2141     Int_t nVertices = nV0s + nCascades + nKinks;
2142
2143     aod->ResetStd(nTracks, nVertices);
2144     AliAODTrack *aodTrack;
2145
2146
2147     // Array to take into account the tracks already added to the AOD
2148     Bool_t * usedTrack = NULL;
2149     if (nTracks>0) {
2150       usedTrack = new Bool_t[nTracks];
2151       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2152     }
2153     // Array to take into account the V0s already added to the AOD
2154     Bool_t * usedV0 = NULL;
2155     if (nV0s>0) {
2156       usedV0 = new Bool_t[nV0s];
2157       for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2158     }
2159     // Array to take into account the kinks already added to the AOD
2160     Bool_t * usedKink = NULL;
2161     if (nKinks>0) {
2162       usedKink = new Bool_t[nKinks];
2163       for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2164     }
2165
2166     // Access to the AOD container of vertices
2167     TClonesArray &vertices = *(aod->GetVertices());
2168     Int_t jVertices=0;
2169
2170     // Access to the AOD container of tracks
2171     TClonesArray &tracks = *(aod->GetTracks());
2172     Int_t jTracks=0;
2173
2174     // Add primary vertex. The primary tracks will be defined
2175     // after the loops on the composite objects (V0, cascades, kinks)
2176     const AliESDVertex *vtx = esd->GetPrimaryVertex();
2177
2178     vtx->GetXYZ(pos); // position
2179     vtx->GetCovMatrix(covVtx); //covariance matrix
2180
2181     AliAODVertex * primary = new(vertices[jVertices++])
2182       AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2183
2184     // Create vertices starting from the most complex objects
2185
2186     // Cascades
2187     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2188       AliESDcascade *cascade = esd->GetCascade(nCascade);
2189
2190       cascade->GetXYZ(pos[0], pos[1], pos[2]);
2191       cascade->GetPosCovXi(covVtx);
2192
2193       // Add the cascade vertex
2194       AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2195                                                                         covVtx,
2196                                                                         cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2197                                                                         primary,
2198                                                                         AliAODVertex::kCascade);
2199
2200       primary->AddDaughter(vcascade);
2201
2202       // Add the V0 from the cascade. The ESD class have to be optimized...
2203       // Now we have to search for the corresponding Vo in the list of V0s
2204       // using the indeces of the positive and negative tracks
2205
2206       Int_t posFromV0 = cascade->GetPindex();
2207       Int_t negFromV0 = cascade->GetNindex();
2208
2209
2210       AliESDv0 * v0 = 0x0;
2211       Int_t indV0 = -1;
2212
2213       for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2214
2215         v0 = esd->GetV0(iV0);
2216         Int_t posV0 = v0->GetPindex();
2217         Int_t negV0 = v0->GetNindex();
2218
2219         if (posV0==posFromV0 && negV0==negFromV0) {
2220           indV0 = iV0;
2221           break;
2222         }
2223       }
2224
2225       AliAODVertex * vV0FromCascade = 0x0;
2226
2227       if (indV0>-1 && !usedV0[indV0] ) {
2228
2229         // the V0 exists in the array of V0s and is not used
2230
2231         usedV0[indV0] = kTRUE;
2232
2233         v0->GetXYZ(pos[0], pos[1], pos[2]);
2234         v0->GetPosCov(covVtx);
2235
2236         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2237                                                                  covVtx,
2238                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2239                                                                  vcascade,
2240                                                                  AliAODVertex::kV0);
2241       } else {
2242
2243         // the V0 doesn't exist in the array of V0s or was used
2244         cerr << "Error: event " << iEvent << " cascade " << nCascade
2245              << " The V0 " << indV0
2246              << " doesn't exist in the array of V0s or was used!" << endl;
2247
2248         cascade->GetXYZ(pos[0], pos[1], pos[2]);
2249         cascade->GetPosCov(covVtx);
2250
2251         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2252                                                                  covVtx,
2253                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2254                                                                  vcascade,
2255                                                                  AliAODVertex::kV0);
2256         vcascade->AddDaughter(vV0FromCascade);
2257       }
2258
2259       // Add the positive tracks from the V0
2260
2261       if (! usedTrack[posFromV0]) {
2262
2263         usedTrack[posFromV0] = kTRUE;
2264
2265         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2266         esdTrack->GetPxPyPz(p);
2267         esdTrack->GetXYZ(pos);
2268         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2269         esdTrack->GetESDpid(pid);
2270
2271         vV0FromCascade->AddDaughter(aodTrack =
2272                                     new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2273                                            esdTrack->GetLabel(),
2274                                            p,
2275                                            kTRUE,
2276                                            pos,
2277                                            kFALSE,
2278                                            covTr,
2279                                            (Short_t)esdTrack->GetSign(),
2280                                            esdTrack->GetITSClusterMap(),
2281                                            pid,
2282                                            vV0FromCascade,
2283                                            kTRUE,  // check if this is right
2284                                            kFALSE, // check if this is right
2285                                            AliAODTrack::kSecondary)
2286                 );
2287         aodTrack->ConvertAliPIDtoAODPID();
2288       }
2289       else {
2290         cerr << "Error: event " << iEvent << " cascade " << nCascade
2291              << " track " << posFromV0 << " has already been used!" << endl;
2292       }
2293
2294       // Add the negative tracks from the V0
2295
2296       if (!usedTrack[negFromV0]) {
2297
2298         usedTrack[negFromV0] = kTRUE;
2299
2300         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2301         esdTrack->GetPxPyPz(p);
2302         esdTrack->GetXYZ(pos);
2303         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2304         esdTrack->GetESDpid(pid);
2305
2306         vV0FromCascade->AddDaughter(aodTrack =
2307                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2308                                            esdTrack->GetLabel(),
2309                                            p,
2310                                            kTRUE,
2311                                            pos,
2312                                            kFALSE,
2313                                            covTr,
2314                                            (Short_t)esdTrack->GetSign(),
2315                                            esdTrack->GetITSClusterMap(),
2316                                            pid,
2317                                            vV0FromCascade,
2318                                            kTRUE,  // check if this is right
2319                                            kFALSE, // check if this is right
2320                                            AliAODTrack::kSecondary)
2321                 );
2322         aodTrack->ConvertAliPIDtoAODPID();
2323       }
2324       else {
2325         cerr << "Error: event " << iEvent << " cascade " << nCascade
2326              << " track " << negFromV0 << " has already been used!" << endl;
2327       }
2328
2329       // Add the bachelor track from the cascade
2330
2331       Int_t bachelor = cascade->GetBindex();
2332
2333       if(!usedTrack[bachelor]) {
2334
2335         usedTrack[bachelor] = kTRUE;
2336
2337         AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2338         esdTrack->GetPxPyPz(p);
2339         esdTrack->GetXYZ(pos);
2340         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2341         esdTrack->GetESDpid(pid);
2342
2343         vcascade->AddDaughter(aodTrack =
2344                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2345                                            esdTrack->GetLabel(),
2346                                            p,
2347                                            kTRUE,
2348                                            pos,
2349                                            kFALSE,
2350                                            covTr,
2351                                            (Short_t)esdTrack->GetSign(),
2352                                            esdTrack->GetITSClusterMap(),
2353                                            pid,
2354                                            vcascade,
2355                                            kTRUE,  // check if this is right
2356                                            kFALSE, // check if this is right
2357                                            AliAODTrack::kSecondary)
2358                 );
2359         aodTrack->ConvertAliPIDtoAODPID();
2360      }
2361       else {
2362         cerr << "Error: event " << iEvent << " cascade " << nCascade
2363              << " track " << bachelor << " has already been used!" << endl;
2364       }
2365
2366       // Add the primary track of the cascade (if any)
2367
2368     } // end of the loop on cascades
2369
2370     // V0s
2371
2372     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2373
2374       if (usedV0[nV0]) continue; // skip if aready added to the AOD
2375
2376       AliESDv0 *v0 = esd->GetV0(nV0);
2377
2378       v0->GetXYZ(pos[0], pos[1], pos[2]);
2379       v0->GetPosCov(covVtx);
2380
2381       AliAODVertex * vV0 =
2382         new(vertices[jVertices++]) AliAODVertex(pos,
2383                                                 covVtx,
2384                                                 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2385                                                 primary,
2386                                                 AliAODVertex::kV0);
2387       primary->AddDaughter(vV0);
2388
2389       Int_t posFromV0 = v0->GetPindex();
2390       Int_t negFromV0 = v0->GetNindex();
2391
2392       // Add the positive tracks from the V0
2393
2394       if (!usedTrack[posFromV0]) {
2395
2396         usedTrack[posFromV0] = kTRUE;
2397
2398         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2399         esdTrack->GetPxPyPz(p);
2400         esdTrack->GetXYZ(pos);
2401         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2402         esdTrack->GetESDpid(pid);
2403
2404         vV0->AddDaughter(aodTrack =
2405                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2406                                            esdTrack->GetLabel(),
2407                                            p,
2408                                            kTRUE,
2409                                            pos,
2410                                            kFALSE,
2411                                            covTr,
2412                                            (Short_t)esdTrack->GetSign(),
2413                                            esdTrack->GetITSClusterMap(),
2414                                            pid,
2415                                            vV0,
2416                                            kTRUE,  // check if this is right
2417                                            kFALSE, // check if this is right
2418                                            AliAODTrack::kSecondary)
2419                 );
2420         aodTrack->ConvertAliPIDtoAODPID();
2421       }
2422       else {
2423         cerr << "Error: event " << iEvent << " V0 " << nV0
2424              << " track " << posFromV0 << " has already been used!" << endl;
2425       }
2426
2427       // Add the negative tracks from the V0
2428
2429       if (!usedTrack[negFromV0]) {
2430
2431         usedTrack[negFromV0] = kTRUE;
2432
2433         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2434         esdTrack->GetPxPyPz(p);
2435         esdTrack->GetXYZ(pos);
2436         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2437         esdTrack->GetESDpid(pid);
2438
2439         vV0->AddDaughter(aodTrack =
2440                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2441                                            esdTrack->GetLabel(),
2442                                            p,
2443                                            kTRUE,
2444                                            pos,
2445                                            kFALSE,
2446                                            covTr,
2447                                            (Short_t)esdTrack->GetSign(),
2448                                            esdTrack->GetITSClusterMap(),
2449                                            pid,
2450                                            vV0,
2451                                            kTRUE,  // check if this is right
2452                                            kFALSE, // check if this is right
2453                                            AliAODTrack::kSecondary)
2454                 );
2455         aodTrack->ConvertAliPIDtoAODPID();
2456       }
2457       else {
2458         cerr << "Error: event " << iEvent << " V0 " << nV0
2459              << " track " << negFromV0 << " has already been used!" << endl;
2460       }
2461
2462     } // end of the loop on V0s
2463
2464     // Kinks: it is a big mess the access to the information in the kinks
2465     // The loop is on the tracks in order to find the mother and daugther of each kink
2466
2467
2468     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2469
2470
2471       AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2472
2473       Int_t ikink = esdTrack->GetKinkIndex(0);
2474
2475       if (ikink) {
2476         // Negative kink index: mother, positive: daughter
2477
2478         // Search for the second track of the kink
2479
2480         for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2481
2482           AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2483
2484           Int_t jkink = esdTrack1->GetKinkIndex(0);
2485
2486           if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2487
2488             // The two tracks are from the same kink
2489
2490             if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2491
2492             Int_t imother = -1;
2493             Int_t idaughter = -1;
2494
2495             if (ikink<0 && jkink>0) {
2496
2497               imother = iTrack;
2498               idaughter = jTrack;
2499             }
2500             else if (ikink>0 && jkink<0) {
2501
2502               imother = jTrack;
2503               idaughter = iTrack;
2504             }
2505             else {
2506               cerr << "Error: Wrong combination of kink indexes: "
2507               << ikink << " " << jkink << endl;
2508               continue;
2509             }
2510
2511             // Add the mother track
2512
2513             AliAODTrack * mother = NULL;
2514
2515             if (!usedTrack[imother]) {
2516
2517               usedTrack[imother] = kTRUE;
2518
2519               AliESDtrack *esdTrack = esd->GetTrack(imother);
2520               esdTrack->GetPxPyPz(p);
2521               esdTrack->GetXYZ(pos);
2522               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2523               esdTrack->GetESDpid(pid);
2524
2525               mother =
2526                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2527                                            esdTrack->GetLabel(),
2528                                            p,
2529                                            kTRUE,
2530                                            pos,
2531                                            kFALSE,
2532                                            covTr,
2533                                            (Short_t)esdTrack->GetSign(),
2534                                            esdTrack->GetITSClusterMap(),
2535                                            pid,
2536                                            primary,
2537                                            kTRUE, // check if this is right
2538                                            kTRUE, // check if this is right
2539                                            AliAODTrack::kPrimary);
2540               primary->AddDaughter(mother);
2541               mother->ConvertAliPIDtoAODPID();
2542             }
2543             else {
2544               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2545               << " track " << imother << " has already been used!" << endl;
2546             }
2547
2548             // Add the kink vertex
2549             AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2550
2551             AliAODVertex * vkink =
2552             new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2553                                                     NULL,
2554                                                     0.,
2555                                                     mother,
2556                                                     AliAODVertex::kKink);
2557             // Add the daughter track
2558
2559             AliAODTrack * daughter = NULL;
2560
2561             if (!usedTrack[idaughter]) {
2562
2563               usedTrack[idaughter] = kTRUE;
2564
2565               AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2566               esdTrack->GetPxPyPz(p);
2567               esdTrack->GetXYZ(pos);
2568               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2569               esdTrack->GetESDpid(pid);
2570
2571               daughter =
2572                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2573                                            esdTrack->GetLabel(),
2574                                            p,
2575                                            kTRUE,
2576                                            pos,
2577                                            kFALSE,
2578                                            covTr,
2579                                            (Short_t)esdTrack->GetSign(),
2580                                            esdTrack->GetITSClusterMap(),
2581                                            pid,
2582                                            vkink,
2583                                            kTRUE, // check if this is right
2584                                            kTRUE, // check if this is right
2585                                            AliAODTrack::kPrimary);
2586               vkink->AddDaughter(daughter);
2587               daughter->ConvertAliPIDtoAODPID();
2588             }
2589             else {
2590               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2591               << " track " << idaughter << " has already been used!" << endl;
2592             }
2593
2594
2595           }
2596         }
2597
2598       }
2599
2600     }
2601
2602
2603     // Tracks (primary and orphan)
2604
2605     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2606
2607
2608       if (usedTrack[nTrack]) continue;
2609
2610       AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2611       esdTrack->GetPxPyPz(p);
2612       esdTrack->GetXYZ(pos);
2613       esdTrack->GetCovarianceXYZPxPyPz(covTr);
2614       esdTrack->GetESDpid(pid);
2615
2616       Float_t impactXY, impactZ;
2617
2618       esdTrack->GetImpactParameters(impactXY,impactZ);
2619
2620       if (impactXY<3) {
2621         // track inside the beam pipe
2622
2623         primary->AddDaughter(aodTrack =
2624             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2625                                          esdTrack->GetLabel(),
2626                                          p,
2627                                          kTRUE,
2628                                          pos,
2629                                          kFALSE,
2630                                          covTr,
2631                                          (Short_t)esdTrack->GetSign(),
2632                                          esdTrack->GetITSClusterMap(),
2633                                          pid,
2634                                          primary,
2635                                          kTRUE, // check if this is right
2636                                          kTRUE, // check if this is right
2637                                          AliAODTrack::kPrimary)
2638             );
2639         aodTrack->ConvertAliPIDtoAODPID();
2640       }
2641       else {
2642         // outside the beam pipe: orphan track
2643             aodTrack =
2644             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2645                                          esdTrack->GetLabel(),
2646                                          p,
2647                                          kTRUE,
2648                                          pos,
2649                                          kFALSE,
2650                                          covTr,
2651                                          (Short_t)esdTrack->GetSign(),
2652                                          esdTrack->GetITSClusterMap(),
2653                                          pid,
2654                                          NULL,
2655                                          kFALSE, // check if this is right
2656                                          kFALSE, // check if this is right
2657                                          AliAODTrack::kOrphan);
2658             aodTrack->ConvertAliPIDtoAODPID();
2659       }
2660     } // end of loop on tracks
2661
2662     // muon tracks
2663     Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2664     for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2665
2666       AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2667       p[0] = esdMuTrack->Px();
2668       p[1] = esdMuTrack->Py();
2669       p[2] = esdMuTrack->Pz();
2670       pos[0] = primary->GetX();
2671       pos[1] = primary->GetY();
2672       pos[2] = primary->GetZ();
2673
2674       // has to be changed once the muon pid is provided by the ESD
2675       for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2676
2677       primary->AddDaughter(
2678           new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2679                                              0, // no label provided
2680                                              p,
2681                                              kTRUE,
2682                                              pos,
2683                                              kFALSE,
2684                                              NULL, // no covariance matrix provided
2685                                              (Short_t)-99, // no charge provided
2686                                              0, // no ITSClusterMap
2687                                              pid,
2688                                              primary,
2689                                              kTRUE,  // check if this is right
2690                                              kTRUE,  // not used for vertex fit
2691                                              AliAODTrack::kPrimary)
2692           );
2693     }
2694
2695     // Access to the AOD container of clusters
2696     TClonesArray &clusters = *(aod->GetClusters());
2697     Int_t jClusters=0;
2698
2699     // Calo Clusters
2700     Int_t nClusters    = esd->GetNumberOfCaloClusters();
2701
2702     for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2703
2704       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2705
2706       Int_t id = cluster->GetID();
2707       Int_t label = -1;
2708       Float_t energy = cluster->GetClusterEnergy();
2709       cluster->GetGlobalPosition(posF);
2710       AliAODVertex *prodVertex = primary;
2711       AliAODTrack *primTrack = NULL;
2712       Char_t ttype=AliAODCluster::kUndef;
2713
2714       if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2715       else if (cluster->IsEMCAL()) {
2716
2717         if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2718           ttype = AliAODCluster::kEMCALPseudoCluster;
2719         else
2720           ttype = AliAODCluster::kEMCALClusterv1;
2721
2722       }
2723
2724       new(clusters[jClusters++]) AliAODCluster(id,
2725                                                label,
2726                                                energy,
2727                                                pos,
2728                                                NULL, // no covariance matrix provided
2729                                                NULL, // no pid for clusters provided
2730                                                prodVertex,
2731                                                primTrack,
2732                                                ttype);
2733
2734     } // end of loop on calo clusters
2735
2736     delete [] usedTrack;
2737     delete [] usedV0;
2738     delete [] usedKink;
2739
2740     // fill the tree for this event
2741     aodTree->Fill();
2742   } // end of event loop
2743
2744   aodTree->GetUserInfo()->Add(aod);
2745
2746   // close ESD file
2747   esdFile->Close();
2748
2749   // write the tree to the specified file
2750   aodFile = aodTree->GetCurrentFile();
2751   aodFile->cd();
2752   aodTree->Write();
2753
2754   return;
2755 }
2756
2757
2758 void AliReconstruction::WriteAlignmentData(AliESD* esd)
2759 {
2760   // Write space-points which are then used in the alignment procedures
2761   // For the moment only ITS, TRD and TPC
2762
2763   // Load TOF clusters
2764   if (fTracker[3]){
2765     fLoader[3]->LoadRecPoints("read");
2766     TTree* tree = fLoader[3]->TreeR();
2767     if (!tree) {
2768       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2769       return;
2770     }
2771     fTracker[3]->LoadClusters(tree);
2772   }
2773   Int_t ntracks = esd->GetNumberOfTracks();
2774   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2775     {
2776       AliESDtrack *track = esd->GetTrack(itrack);
2777       Int_t nsp = 0;
2778       Int_t idx[200];
2779       for (Int_t iDet = 3; iDet >= 0; iDet--)
2780         nsp += track->GetNcls(iDet);
2781       if (nsp) {
2782         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2783         track->SetTrackPointArray(sp);
2784         Int_t isptrack = 0;
2785         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2786           AliTracker *tracker = fTracker[iDet];
2787           if (!tracker) continue;
2788           Int_t nspdet = track->GetNcls(iDet);
2789           if (nspdet <= 0) continue;
2790           track->GetClusters(iDet,idx);
2791           AliTrackPoint p;
2792           Int_t isp = 0;
2793           Int_t isp2 = 0;
2794           while (isp < nspdet) {
2795             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2796             const Int_t kNTPCmax = 159;
2797             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2798             if (!isvalid) continue;
2799             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2800           }
2801         }
2802       }
2803     }
2804   if (fTracker[3]){
2805     fTracker[3]->UnloadClusters();
2806     fLoader[3]->UnloadRecPoints();
2807   }
2808 }
2809
2810 //_____________________________________________________________________________
2811 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2812 {
2813   // The method reads the raw-data error log
2814   // accumulated within the rawReader.
2815   // It extracts the raw-data errors related to
2816   // the current event and stores them into
2817   // a TClonesArray inside the esd object.
2818
2819   if (!fRawReader) return;
2820
2821   for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2822
2823     AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2824     if (!log) continue;
2825     if (iEvent != log->GetEventNumber()) continue;
2826
2827     esd->AddRawDataErrorLog(log);
2828   }
2829
2830 }