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