STEER/AliReconstruction.cxx

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