STEER/AliReconstruction.cxx

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