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