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