STEER/AliReconstruction.cxx

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