X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=STEER%2FAliReconstruction.cxx;h=a9ee0e1073a003a5ceaca136d2b651b24f4698cb;hp=ea359524f2a98d4c45979cd93d3ee2a22010ba29;hb=1f46a9aec2e119240d1d2960672862a073fba4da;hpb=85555ca82828c65e6569c303fc08111f25f333b6 diff --git a/STEER/AliReconstruction.cxx b/STEER/AliReconstruction.cxx index ea359524f2a..a9ee0e1073a 100644 --- a/STEER/AliReconstruction.cxx +++ b/STEER/AliReconstruction.cxx @@ -26,74 +26,177 @@ // rec.Run(); // // // // The Run method returns kTRUE in case of successful execution. // +// // +// If the input to the reconstruction are not simulated digits but raw data, // +// this can be specified by an argument of the Run method or by the method // +// // +// rec.SetInput("..."); // +// // +// The input formats and the corresponding argument are: // +// - DDL raw data files: directory name, ends with "/" // +// - raw data root file: root file name, extension ".root" // +// - raw data DATE file: DATE file name, any other non-empty string // +// - MC root files : empty string, default // +// // +// By default all events are reconstructed. The reconstruction can be // +// limited to a range of events by giving the index of the first and the // +// last event as an argument to the Run method or by calling // +// // +// rec.SetEventRange(..., ...); // +// // +// The index -1 (default) can be used for the last event to indicate no // +// upper limit of the event range. // +// // // The name of the galice file can be changed from the default // -// "galice.root" by // +// "galice.root" by passing it as argument to the AliReconstruction // +// constructor or by // // // // rec.SetGAliceFile("..."); // // // -// The reconstruction can be switched on or off for individual detectors by // +// The local reconstruction can be switched on or off for individual // +// detectors by // // // -// rec.SetRunReconstruction("..."); // +// rec.SetRunLocalReconstruction("..."); // // // // The argument is a (case sensitive) string with the names of the // // detectors separated by a space. The special string "ALL" selects all // // available detectors. This is the default. // // // -// The tracking in ITS, TPC and TRD and the creation of ESD tracks can be // -// switched off by // +// The reconstruction of the primary vertex position can be switched off by // +// // +// rec.SetRunVertexFinder(kFALSE); // // // -// rec.SetRunTracking(kFALSE); // +// The tracking and the creation of ESD tracks can be switched on for // +// selected detectors by // +// // +// rec.SetRunTracking("..."); // // // // The filling of additional ESD information can be steered by // // // // rec.SetFillESD("..."); // // // -// Again, the string specifies the list of detectors. The default is "ALL". // +// Again, for both methods the string specifies the list of detectors. // +// The default is "ALL". // +// // +// The call of the shortcut method // +// // +// rec.SetRunReconstruction("..."); // +// // +// is equivalent to calling SetRunLocalReconstruction, SetRunTracking and // +// SetFillESD with the same detector selecting string as argument. // +// // +// The reconstruction requires digits or raw data as input. For the creation // +// of digits and raw data have a look at the class AliSimulation. // // // -// The reconstruction requires digits as input. For the creation of digits // -// have a look at the class AliSimulation. // +// For debug purposes the method SetCheckPointLevel can be used. If the // +// argument is greater than 0, files with ESD events will be written after // +// selected steps of the reconstruction for each event: // +// level 1: after tracking and after filling of ESD (final) // +// level 2: in addition after each tracking step // +// level 3: in addition after the filling of ESD for each detector // +// If a final check point file exists for an event, this event will be // +// skipped in the reconstruction. The tracking and the filling of ESD for // +// a detector will be skipped as well, if the corresponding check point // +// file exists. The ESD event will then be loaded from the file instead. // // // /////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include +#include +#include #include "AliReconstruction.h" +#include "AliReconstructor.h" +#include "AliLog.h" #include "AliRunLoader.h" #include "AliRun.h" -#include "AliModule.h" -#include "AliDetector.h" +#include "AliRawReaderFile.h" +#include "AliRawReaderDate.h" +#include "AliRawReaderRoot.h" #include "AliTracker.h" #include "AliESD.h" +#include "AliESDVertex.h" +#include "AliVertexer.h" #include "AliHeader.h" #include "AliGenEventHeader.h" +#include "AliPID.h" #include "AliESDpid.h" -#include - +#include "AliMagF.h" ClassImp(AliReconstruction) //_____________________________________________________________________________ -AliReconstruction::AliReconstruction(const char* name, const char* title) : - TNamed(name, title) +const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"}; + +//_____________________________________________________________________________ +AliReconstruction::AliReconstruction(const char* gAliceFilename, + const char* name, const char* title) : + TNamed(name, title), + + fRunLocalReconstruction("ALL"), + fRunVertexFinder(kTRUE), + fRunHLTTracking(kFALSE), + fRunTracking("ALL"), + fFillESD("ALL"), + fGAliceFileName(gAliceFilename), + fInput(""), + fFirstEvent(0), + fLastEvent(-1), + fStopOnError(kFALSE), + fCheckPointLevel(0), + fOptions(), + + fRunLoader(NULL), + fRawReader(NULL), + + fVertexer(NULL) { // create reconstruction object with default parameters - - Init(); + + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + fReconstructor[iDet] = NULL; + fLoader[iDet] = NULL; + fTracker[iDet] = NULL; + } + // AliPID::Init(); } //_____________________________________________________________________________ AliReconstruction::AliReconstruction(const AliReconstruction& rec) : - TNamed(rec) + TNamed(rec), + + fRunLocalReconstruction(rec.fRunLocalReconstruction), + fRunVertexFinder(rec.fRunVertexFinder), + fRunHLTTracking(rec.fRunHLTTracking), + fRunTracking(rec.fRunTracking), + fFillESD(rec.fFillESD), + fGAliceFileName(rec.fGAliceFileName), + fInput(rec.fInput), + fFirstEvent(rec.fFirstEvent), + fLastEvent(rec.fLastEvent), + fStopOnError(rec.fStopOnError), + fCheckPointLevel(0), + fOptions(), + + fRunLoader(NULL), + fRawReader(NULL), + + fVertexer(NULL) { // copy constructor - fRunReconstruction = rec.fRunReconstruction; - fRunTracking = rec.fRunTracking; - fStopOnError = rec.fStopOnError; - - fGAliceFileName = rec.fGAliceFileName; - - fRunLoader = NULL; + for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) { + if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone()); + } + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + fReconstructor[iDet] = NULL; + fLoader[iDet] = NULL; + fTracker[iDet] = NULL; + } } //_____________________________________________________________________________ @@ -111,362 +214,609 @@ AliReconstruction::~AliReconstruction() { // clean up + CleanUp(); + fOptions.Delete(); } + //_____________________________________________________________________________ -void AliReconstruction::Init() +void AliReconstruction::SetGAliceFile(const char* fileName) { -// set default parameters - - fRunReconstruction = "ALL"; - fRunTracking = kTRUE; - fFillESD = "ALL"; - fStopOnError = kFALSE; - - fGAliceFileName = "galice.root"; +// set the name of the galice file - fRunLoader = NULL; + fGAliceFileName = fileName; } - //_____________________________________________________________________________ -void AliReconstruction::SetGAliceFile(const char* fileName) +void AliReconstruction::SetOption(const char* detector, const char* option) { -// set the name of the galice file +// set options for the reconstruction of a detector - fGAliceFileName = fileName; + TObject* obj = fOptions.FindObject(detector); + if (obj) fOptions.Remove(obj); + fOptions.Add(new TNamed(detector, option)); } //_____________________________________________________________________________ -Bool_t AliReconstruction::Run() +Bool_t AliReconstruction::Run(const char* input, + Int_t firstEvent, Int_t lastEvent) { // run the reconstruction - // open the run loader - if (fRunLoader) delete fRunLoader; - fRunLoader = AliRunLoader::Open(fGAliceFileName.Data()); - if (!fRunLoader) { - Error("Run", "no run loader found in file %s", - fGAliceFileName.Data()); - return kFALSE; - } - fRunLoader->LoadgAlice(); - gAlice = fRunLoader->GetAliRun(); - if (!gAlice) { - Error("Run", "no gAlice object found in file %s", - fGAliceFileName.Data()); - return kFALSE; + // set the input + if (!input) input = fInput.Data(); + TString fileName(input); + if (fileName.EndsWith("/")) { + fRawReader = new AliRawReaderFile(fileName); + } else if (fileName.EndsWith(".root")) { + fRawReader = new AliRawReaderRoot(fileName); + } else if (!fileName.IsNull()) { + fRawReader = new AliRawReaderDate(fileName); + fRawReader->SelectEvents(7); } + // get the run loader + if (!InitRunLoader()) return kFALSE; + // local reconstruction - if (!fRunReconstruction.IsNull()) { - if (!RunReconstruction(fRunReconstruction)) { - if (fStopOnError) return kFALSE; + if (!fRunLocalReconstruction.IsNull()) { + if (!RunLocalReconstruction(fRunLocalReconstruction)) { + if (fStopOnError) {CleanUp(); return kFALSE;} } } - if (!fRunTracking && fFillESD.IsNull()) return kTRUE; - - // get loaders and trackers - fITSLoader = fRunLoader->GetLoader("ITSLoader"); - if (!fITSLoader) { - Error("Run", "no ITS loader found"); - if (fStopOnError) return kFALSE; - } - fITSTracker = NULL; - if (gAlice->GetDetector("ITS")) { - fITSTracker = gAlice->GetDetector("ITS")->CreateTracker(); - } - if (!fITSTracker) { - Error("Run", "couldn't create a tracker for ITS"); - if (fStopOnError) return kFALSE; - } - - fTPCLoader = fRunLoader->GetLoader("TPCLoader"); - if (!fTPCLoader) { - Error("Run", "no TPC loader found"); - if (fStopOnError) return kFALSE; - } - fTPCTracker = NULL; - if (gAlice->GetDetector("TPC")) { - fTPCTracker = gAlice->GetDetector("TPC")->CreateTracker(); - } - if (!fTPCTracker) { - Error("Run", "couldn't create a tracker for TPC"); - if (fStopOnError) return kFALSE; +// if (!fRunVertexFinder && fRunTracking.IsNull() && +// fFillESD.IsNull()) return kTRUE; + + // get vertexer + if (fRunVertexFinder && !CreateVertexer()) { + if (fStopOnError) { + CleanUp(); + return kFALSE; + } } - fTRDLoader = fRunLoader->GetLoader("TRDLoader"); - if (!fTRDLoader) { - Error("Run", "no TRD loader found"); - if (fStopOnError) return kFALSE; - } - fTRDTracker = NULL; - if (gAlice->GetDetector("TRD")) { - fTRDTracker = gAlice->GetDetector("TRD")->CreateTracker(); - } - if (!fTRDTracker) { - Error("Run", "couldn't create a tracker for TRD"); - if (fStopOnError) return kFALSE; + // get trackers + if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) { + if (fStopOnError) { + CleanUp(); + return kFALSE; + } } - fTOFLoader = fRunLoader->GetLoader("TOFLoader"); - if (!fTOFLoader) { - Error("Run", "no TOF loader found"); - if (fStopOnError) return kFALSE; - } - fTOFTracker = NULL; - if (gAlice->GetDetector("TOF")) { - fTOFTracker = gAlice->GetDetector("TOF")->CreateTracker(); - } - if (!fTOFTracker) { - Error("Run", "couldn't create a tracker for TOF"); - if (fStopOnError) return kFALSE; + // get the possibly already existing ESD file and tree + AliESD* esd = new AliESD; AliESD* hltesd = new AliESD; + TFile* fileOld = NULL; + TTree* treeOld = NULL; TTree *hlttreeOld = NULL; + if (!gSystem->AccessPathName("AliESDs.root")){ + gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE); + fileOld = TFile::Open("AliESDs.old.root"); + if (fileOld && fileOld->IsOpen()) { + treeOld = (TTree*) fileOld->Get("esdTree"); + if (treeOld) treeOld->SetBranchAddress("ESD", &esd); + hlttreeOld = (TTree*) fileOld->Get("HLTesdTree"); + if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd); + } } - // create the ESD output file + // create the ESD output file and tree TFile* file = TFile::Open("AliESDs.root", "RECREATE"); if (!file->IsOpen()) { - Error("Run", "opening AliESDs.root failed"); - if (fStopOnError) return kFALSE; + AliError("opening AliESDs.root failed"); + if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;} } + TTree* tree = new TTree("esdTree", "Tree with ESD objects"); + tree->Branch("ESD", "AliESD", &esd); + TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects"); + hlttree->Branch("ESD", "AliESD", &hltesd); + delete esd; delete hltesd; + esd = NULL; hltesd = NULL; + gROOT->cd(); // loop over events + if (fRawReader) fRawReader->RewindEvents(); + for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) { - Info("Run", "processing event %d", iEvent); - AliESD* esd = new AliESD; + if (fRawReader) fRawReader->NextEvent(); + if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) { + // copy old ESD to the new one + if (treeOld) { + treeOld->SetBranchAddress("ESD", &esd); + treeOld->GetEntry(iEvent); + } + tree->Fill(); + if (hlttreeOld) { + hlttreeOld->SetBranchAddress("ESD", &hltesd); + hlttreeOld->GetEntry(iEvent); + } + hlttree->Fill(); + continue; + } + + AliInfo(Form("processing event %d", iEvent)); fRunLoader->GetEvent(iEvent); - esd->SetRunNumber(gAlice->GetRunNumber()); - esd->SetEventNumber(gAlice->GetEvNumber()); + + char fileName[256]; + sprintf(fileName, "ESD_%d.%d_final.root", + fRunLoader->GetHeader()->GetRun(), + fRunLoader->GetHeader()->GetEventNrInRun()); + if (!gSystem->AccessPathName(fileName)) continue; + + // local reconstruction + if (!fRunLocalReconstruction.IsNull()) { + if (!RunLocalEventReconstruction(fRunLocalReconstruction)) { + if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;} + } + } + + esd = new AliESD; hltesd = new AliESD; + esd->SetRunNumber(fRunLoader->GetHeader()->GetRun()); + hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun()); + esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun()); + hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun()); + if (gAlice) { + esd->SetMagneticField(gAlice->Field()->SolenoidField()); + hltesd->SetMagneticField(gAlice->Field()->SolenoidField()); + } else { + // ??? + } + + // vertex finder + if (fRunVertexFinder) { + if (!ReadESD(esd, "vertex")) { + if (!RunVertexFinder(esd)) { + if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;} + } + if (fCheckPointLevel > 0) WriteESD(esd, "vertex"); + } + } + + // HLT tracking + if (!fRunTracking.IsNull()) { + if (fRunHLTTracking) { + hltesd->SetVertex(esd->GetVertex()); + if (!RunHLTTracking(hltesd)) { + if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;} + } + } + } // barrel tracking - if (fRunTracking) { - if (!RunTracking(esd)) { - if (fStopOnError) return kFALSE; + if (!fRunTracking.IsNull()) { + if (!ReadESD(esd, "tracking")) { + if (!RunTracking(esd)) { + if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;} + } + if (fCheckPointLevel > 0) WriteESD(esd, "tracking"); } } // fill ESD if (!fFillESD.IsNull()) { if (!FillESD(esd, fFillESD)) { - if (fStopOnError) return kFALSE; + if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;} } } // combined PID AliESDpid::MakePID(esd); + if (fCheckPointLevel > 1) WriteESD(esd, "PID"); // write ESD - char name[100]; - sprintf(name, "ESD%d", iEvent); - file->cd(); - if (!esd->Write(name)) { - Error("Run", "writing ESD failed"); - if (fStopOnError) return kFALSE; - } + tree->Fill(); + // write HLT ESD + hlttree->Fill(); + + if (fCheckPointLevel > 0) WriteESD(esd, "final"); + delete esd; delete hltesd; + esd = NULL; hltesd = NULL; } - file->Close(); + file->cd(); + tree->Write(); + hlttree->Write(); + CleanUp(file, fileOld); return kTRUE; } //_____________________________________________________________________________ -Bool_t AliReconstruction::RunReconstruction(const TString& detectors) +Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors) { -// run the reconstruction +// run the local reconstruction + + TStopwatch stopwatch; + stopwatch.Start(); + + TString detStr = detectors; + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + if (!IsSelected(fgkDetectorName[iDet], detStr)) continue; + AliReconstructor* reconstructor = GetReconstructor(iDet); + if (!reconstructor) continue; + if (reconstructor->HasLocalReconstruction()) continue; + + AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet])); + TStopwatch stopwatchDet; + stopwatchDet.Start(); + if (fRawReader) { + fRawReader->RewindEvents(); + reconstructor->Reconstruct(fRunLoader, fRawReader); + } else { + reconstructor->Reconstruct(fRunLoader); + } + AliInfo(Form("execution time for %s:", fgkDetectorName[iDet])); + ToAliInfo(stopwatchDet.Print()); + } + + if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) { + AliError(Form("the following detectors were not found: %s", + detStr.Data())); + if (fStopOnError) return kFALSE; + } + + AliInfo("execution time:"); + ToAliInfo(stopwatch.Print()); + + return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors) +{ +// run the local reconstruction TStopwatch stopwatch; stopwatch.Start(); TString detStr = detectors; - TObjArray* detArray = gAlice->Detectors(); - for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) { - AliModule* det = (AliModule*) detArray->At(iDet); - if (!det || !det->IsActive()) continue; - if (IsSelected(det->GetName(), detStr)) { - Info("RunReconstruction", "running reconstruction for %s", - det->GetName()); + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + if (!IsSelected(fgkDetectorName[iDet], detStr)) continue; + AliReconstructor* reconstructor = GetReconstructor(iDet); + if (!reconstructor) continue; + AliLoader* loader = fLoader[iDet]; + + // conversion of digits + if (fRawReader && reconstructor->HasDigitConversion()) { + AliInfo(Form("converting raw data digits into root objects for %s", + fgkDetectorName[iDet])); TStopwatch stopwatchDet; stopwatchDet.Start(); - det->Reconstruct(); - Info("RunReconstruction", "execution time for %s:", det->GetName()); - stopwatchDet.Print(); + loader->LoadDigits("update"); + loader->CleanDigits(); + loader->MakeDigitsContainer(); + TTree* digitsTree = loader->TreeD(); + reconstructor->ConvertDigits(fRawReader, digitsTree); + loader->WriteDigits("OVERWRITE"); + loader->UnloadDigits(); + AliDebug(1, Form("execution time for %s:", fgkDetectorName[iDet])); + ToAliDebug(1, stopwatchDet.Print()); } + + // local reconstruction + if (!reconstructor->HasLocalReconstruction()) continue; + AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet])); + TStopwatch stopwatchDet; + stopwatchDet.Start(); + loader->LoadRecPoints("update"); + loader->CleanRecPoints(); + loader->MakeRecPointsContainer(); + TTree* clustersTree = loader->TreeR(); + if (fRawReader && !reconstructor->HasDigitConversion()) { + reconstructor->Reconstruct(fRawReader, clustersTree); + } else { + loader->LoadDigits("read"); + TTree* digitsTree = loader->TreeD(); + if (!digitsTree) { + AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet])); + if (fStopOnError) return kFALSE; + } else { + reconstructor->Reconstruct(digitsTree, clustersTree); + } + loader->UnloadDigits(); + } + loader->WriteRecPoints("OVERWRITE"); + loader->UnloadRecPoints(); + AliDebug(1, Form("execution time for %s:", fgkDetectorName[iDet])); + ToAliDebug(1, stopwatchDet.Print()); } if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) { - Error("RunReconstruction", "the following detectors were not found: %s", - detStr.Data()); + AliError(Form("the following detectors were not found: %s", + detStr.Data())); if (fStopOnError) return kFALSE; } - Info("RunReconstruction", "execution time:"); - stopwatch.Print(); + AliInfo("execution time:"); + ToAliInfo(stopwatch.Print()); return kTRUE; } //_____________________________________________________________________________ -Bool_t AliReconstruction::RunTracking(AliESD* esd) +Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd) { // run the barrel tracking TStopwatch stopwatch; stopwatch.Start(); - // get the primary vertex (from MC for the moment) - TArrayF vertex(3); - fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex); - Double_t vtxPos[3] = {vertex[0], vertex[1], vertex[2]}; - Double_t vtxCov[6] = { - 0.005, - 0.000, 0.005, - 0.000, 0.000, 0.010 - }; - Double_t vtxErr[3] = {vtxCov[0], vtxCov[2], vtxCov[5]}; // diag. elements - esd->SetVertex(vtxPos, vtxCov); - fITSTracker->SetVertex(vtxPos, vtxErr); - fTPCTracker->SetVertex(vtxPos, vtxErr); - fTRDTracker->SetVertex(vtxPos, vtxErr); - - // TPC tracking - Info("RunTracking", "TPC tracking"); - fTPCLoader->LoadRecPoints("read"); - TTree* tpcTree = fTPCLoader->TreeR(); - if (!tpcTree) { - Error("RunTracking", "Can't get the TPC cluster tree"); - return kFALSE; - } - fTPCTracker->LoadClusters(tpcTree); - if (fTPCTracker->Clusters2Tracks(esd) != 0) { - Error("RunTracking", "TPC Clusters2Tracks failed"); - return kFALSE; + AliESDVertex* vertex = NULL; + Double_t vtxPos[3] = {0, 0, 0}; + Double_t vtxErr[3] = {0.07, 0.07, 0.1}; + TArrayF mcVertex(3); + if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) { + fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex); + for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i]; } - gAlice->GetDetector("TPC")->FillESD(esd); // preliminary PID - AliESDpid::MakePID(esd); // for the ITS tracker + if (fVertexer) { + AliInfo("running the ITS vertex finder"); + if (fLoader[0]) fLoader[0]->LoadRecPoints(); + vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber()); + if (fLoader[0]) fLoader[0]->UnloadRecPoints(); + if(!vertex){ + AliWarning("Vertex not found"); + vertex = new AliESDVertex(); + } + else { + vertex->SetTruePos(vtxPos); // store also the vertex from MC + } - // ITS tracking - Info("RunTracking", "ITS tracking"); - fITSLoader->LoadRecPoints("read"); - TTree* itsTree = fITSLoader->TreeR(); - if (!itsTree) { - Error("RunTracking", "Can't get the ITS cluster tree"); - return kFALSE; - } - fITSTracker->LoadClusters(itsTree); - if (fITSTracker->Clusters2Tracks(esd) != 0) { - Error("RunTracking", "ITS Clusters2Tracks failed"); - return kFALSE; + } else { + AliInfo("getting the primary vertex from MC"); + vertex = new AliESDVertex(vtxPos, vtxErr); } - // ITS back propagation - Info("RunTracking", "ITS back propagation"); - if (fITSTracker->PropagateBack(esd) != 0) { - Error("RunTracking", "ITS backward propagation failed"); - return kFALSE; + if (vertex) { + vertex->GetXYZ(vtxPos); + vertex->GetSigmaXYZ(vtxErr); + } else { + AliWarning("no vertex reconstructed"); + vertex = new AliESDVertex(vtxPos, vtxErr); } + esd->SetVertex(vertex); + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr); + } + delete vertex; - // TPC back propagation - Info("RunTracking", "TPC back propagation"); - if (fTPCTracker->PropagateBack(esd) != 0) { - Error("RunTracking", "TPC backward propagation failed"); - return kFALSE; - } + AliInfo("execution time:"); + ToAliInfo(stopwatch.Print()); - // TRD back propagation - Info("RunTracking", "TRD back propagation"); - fTRDLoader->LoadRecPoints("read"); - TTree* trdTree = fTRDLoader->TreeR(); - if (!trdTree) { - Error("RunTracking", "Can't get the TRD cluster tree"); - return kFALSE; - } - fTRDTracker->LoadClusters(trdTree); - if (fTRDTracker->PropagateBack(esd) != 0) { - Error("RunTracking", "TRD backward propagation failed"); + return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd) +{ +// run the HLT barrel tracking + + TStopwatch stopwatch; + stopwatch.Start(); + + if (!fRunLoader) { + AliError("Missing runLoader!"); return kFALSE; } - // TOF back propagation - Info("RunTracking", "TOF back propagation"); - fTOFLoader->LoadDigits("read"); - TTree* tofTree = fTOFLoader->TreeD(); - if (!tofTree) { - Error("RunTracking", "Can't get the TOF digits tree"); - return kFALSE; - } - fTOFTracker->LoadClusters(tofTree); - if (fTOFTracker->PropagateBack(esd) != 0) { - Error("RunTracking", "TOF backward propagation failed"); - return kFALSE; + AliInfo("running HLT tracking"); + + // Get a pointer to the HLT reconstructor + AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1); + if (!reconstructor) return kFALSE; + + // TPC + ITS + for (Int_t iDet = 1; iDet >= 0; iDet--) { + TString detName = fgkDetectorName[iDet]; + AliDebug(1, Form("%s HLT tracking", detName.Data())); + reconstructor->SetOption(detName.Data()); + AliTracker *tracker = reconstructor->CreateTracker(fRunLoader); + if (!tracker) { + AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data())); + if (fStopOnError) return kFALSE; + } + Double_t vtxPos[3]; + Double_t vtxErr[3]={0.005,0.005,0.010}; + const AliESDVertex *vertex = esd->GetVertex(); + vertex->GetXYZ(vtxPos); + tracker->SetVertex(vtxPos,vtxErr); + if(iDet != 1) { + fLoader[iDet]->LoadRecPoints("read"); + TTree* tree = fLoader[iDet]->TreeR(); + if (!tree) { + AliError(Form("Can't get the %s cluster tree", detName.Data())); + return kFALSE; + } + tracker->LoadClusters(tree); + } + if (tracker->Clusters2Tracks(esd) != 0) { + AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet])); + return kFALSE; + } + if(iDet != 1) { + tracker->UnloadClusters(); + } + delete tracker; } - fTOFTracker->UnloadClusters(); - fTOFLoader->UnloadDigits(); - // TRD inward refit - Info("RunTracking", "TRD inward refit"); - if (fTRDTracker->RefitInward(esd) != 0) { - Error("RunTracking", "TRD inward refit failed"); - return kFALSE; + AliInfo("execution time:"); + ToAliInfo(stopwatch.Print()); + + return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliReconstruction::RunTracking(AliESD*& esd) +{ +// run the barrel tracking + + TStopwatch stopwatch; + stopwatch.Start(); + + AliInfo("running tracking"); + + // pass 1: TPC + ITS inwards + for (Int_t iDet = 1; iDet >= 0; iDet--) { + if (!fTracker[iDet]) continue; + AliDebug(1, Form("%s tracking", fgkDetectorName[iDet])); + + // load clusters + fLoader[iDet]->LoadRecPoints("read"); + TTree* tree = fLoader[iDet]->TreeR(); + if (!tree) { + AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet])); + return kFALSE; + } + fTracker[iDet]->LoadClusters(tree); + + // run tracking + if (fTracker[iDet]->Clusters2Tracks(esd) != 0) { + AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet])); + return kFALSE; + } + if (fCheckPointLevel > 1) { + WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet])); + } + // preliminary PID in TPC needed by the ITS tracker + if (iDet == 1) { + GetReconstructor(1)->FillESD(fRunLoader, esd); + GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd); + AliESDpid::MakePID(esd); + } } - fTRDTracker->UnloadClusters(); - fTRDLoader->UnloadRecPoints(); - - // TPC inward refit - Info("RunTracking", "TPC inward refit"); - if (fTPCTracker->RefitInward(esd) != 0) { - Error("RunTracking", "TPC inward refit failed"); - return kFALSE; + + // pass 2: ALL backwards + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + if (!fTracker[iDet]) continue; + AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet])); + + // load clusters + if (iDet > 1) { // all except ITS, TPC + TTree* tree = NULL; + if (iDet == 3) { // TOF + fLoader[iDet]->LoadDigits("read"); + tree = fLoader[iDet]->TreeD(); + } else { + fLoader[iDet]->LoadRecPoints("read"); + tree = fLoader[iDet]->TreeR(); + } + if (!tree) { + AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet])); + return kFALSE; + } + fTracker[iDet]->LoadClusters(tree); + } + + // run tracking + if (fTracker[iDet]->PropagateBack(esd) != 0) { + AliError(Form("%s backward propagation failed", fgkDetectorName[iDet])); + return kFALSE; + } + if (fCheckPointLevel > 1) { + WriteESD(esd, Form("%s.back", fgkDetectorName[iDet])); + } + + // unload clusters + if (iDet > 2) { // all except ITS, TPC, TRD + fTracker[iDet]->UnloadClusters(); + if (iDet == 3) { // TOF + fLoader[iDet]->UnloadDigits(); + } else { + fLoader[iDet]->UnloadRecPoints(); + } + } } - fTPCTracker->UnloadClusters(); - fTPCLoader->UnloadRecPoints(); - - // ITS inward refit - Info("RunTracking", "ITS inward refit"); - if (fITSTracker->RefitInward(esd) != 0) { - Error("RunTracking", "ITS inward refit failed"); - return kFALSE; + + // pass 3: TRD + TPC + ITS refit inwards + for (Int_t iDet = 2; iDet >= 0; iDet--) { + if (!fTracker[iDet]) continue; + AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet])); + + // run tracking + if (fTracker[iDet]->RefitInward(esd) != 0) { + AliError(Form("%s inward refit failed", fgkDetectorName[iDet])); + return kFALSE; + } + if (fCheckPointLevel > 1) { + WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet])); + } + + // unload clusters + fTracker[iDet]->UnloadClusters(); + fLoader[iDet]->UnloadRecPoints(); } - fITSTracker->UnloadClusters(); - fITSLoader->UnloadRecPoints(); - Info("RunTracking", "execution time:"); - stopwatch.Print(); + AliInfo("execution time:"); + ToAliInfo(stopwatch.Print()); return kTRUE; } //_____________________________________________________________________________ -Bool_t AliReconstruction::FillESD(AliESD* esd, const TString& detectors) +Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors) { // fill the event summary data TStopwatch stopwatch; stopwatch.Start(); + AliInfo("filling ESD"); TString detStr = detectors; - TObjArray* detArray = gAlice->Detectors(); - for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) { - AliModule* det = (AliModule*) detArray->At(iDet); - if (!det || !det->IsActive()) continue; - if (IsSelected(det->GetName(), detStr)) { - Info("FillESD", "filling ESD for %s", - det->GetName()); - det->FillESD(esd); + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + if (!IsSelected(fgkDetectorName[iDet], detStr)) continue; + AliReconstructor* reconstructor = GetReconstructor(iDet); + if (!reconstructor) continue; + + if (!ReadESD(esd, fgkDetectorName[iDet])) { + AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet])); + TTree* clustersTree = NULL; + if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) { + fLoader[iDet]->LoadRecPoints("read"); + clustersTree = fLoader[iDet]->TreeR(); + if (!clustersTree) { + AliError(Form("Can't get the %s clusters tree", + fgkDetectorName[iDet])); + if (fStopOnError) return kFALSE; + } + } + if (fRawReader && !reconstructor->HasDigitConversion()) { + reconstructor->FillESD(fRawReader, clustersTree, esd); + } else { + TTree* digitsTree = NULL; + if (fLoader[iDet]) { + fLoader[iDet]->LoadDigits("read"); + digitsTree = fLoader[iDet]->TreeD(); + if (!digitsTree) { + AliError(Form("Can't get the %s digits tree", + fgkDetectorName[iDet])); + if (fStopOnError) return kFALSE; + } + } + reconstructor->FillESD(digitsTree, clustersTree, esd); + if (fLoader[iDet]) fLoader[iDet]->UnloadDigits(); + } + if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) { + fLoader[iDet]->UnloadRecPoints(); + } + + if (fRawReader) { + reconstructor->FillESD(fRunLoader, fRawReader, esd); + } else { + reconstructor->FillESD(fRunLoader, esd); + } + if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]); } } if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) { - Error("FillESD", "the following detectors were not found: %s", - detStr.Data()); + AliError(Form("the following detectors were not found: %s", + detStr.Data())); if (fStopOnError) return kFALSE; } - Info("FillESD", "execution time:"); - stopwatch.Print(); + AliInfo("execution time:"); + ToAliInfo(stopwatch.Print()); return kTRUE; } @@ -504,3 +854,292 @@ Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const return result; } + +//_____________________________________________________________________________ +Bool_t AliReconstruction::InitRunLoader() +{ +// get or create the run loader + + if (gAlice) delete gAlice; + gAlice = NULL; + + if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists + // load all base libraries to get the loader classes + TString libs = gSystem->GetLibraries(); + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + TString detName = fgkDetectorName[iDet]; + if (detName == "HLT") continue; + if (libs.Contains("lib" + detName + "base.so")) continue; + gSystem->Load("lib" + detName + "base.so"); + } + fRunLoader = AliRunLoader::Open(fGAliceFileName.Data()); + if (!fRunLoader) { + AliError(Form("no run loader found in file %s", fGAliceFileName.Data())); + CleanUp(); + return kFALSE; + } + fRunLoader->CdGAFile(); + if (gFile->GetKey(AliRunLoader::GetGAliceName())) { + if (fRunLoader->LoadgAlice() == 0) { + gAlice = fRunLoader->GetAliRun(); + AliTracker::SetFieldMap(gAlice->Field()); + } + } + if (!gAlice && !fRawReader) { + AliError(Form("no gAlice object found in file %s", + fGAliceFileName.Data())); + CleanUp(); + return kFALSE; + } + + } else { // galice.root does not exist + if (!fRawReader) { + AliError(Form("the file %s does not exist", fGAliceFileName.Data())); + CleanUp(); + return kFALSE; + } + fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(), + AliConfig::GetDefaultEventFolderName(), + "recreate"); + if (!fRunLoader) { + AliError(Form("could not create run loader in file %s", + fGAliceFileName.Data())); + CleanUp(); + return kFALSE; + } + fRunLoader->MakeTree("E"); + Int_t iEvent = 0; + while (fRawReader->NextEvent()) { + fRunLoader->SetEventNumber(iEvent); + fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), + iEvent, iEvent); + fRunLoader->MakeTree("H"); + fRunLoader->TreeE()->Fill(); + iEvent++; + } + fRawReader->RewindEvents(); + fRunLoader->WriteHeader("OVERWRITE"); + fRunLoader->CdGAFile(); + fRunLoader->Write(0, TObject::kOverwrite); +// AliTracker::SetFieldMap(???); + } + + return kTRUE; +} + +//_____________________________________________________________________________ +AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet) +{ +// get the reconstructor object and the loader for a detector + + if (fReconstructor[iDet]) return fReconstructor[iDet]; + + // load the reconstructor object + TPluginManager* pluginManager = gROOT->GetPluginManager(); + TString detName = fgkDetectorName[iDet]; + TString recName = "Ali" + detName + "Reconstructor"; + if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL; + + if (detName == "HLT") { + if (!gROOT->GetClass("AliLevel3")) { + gSystem->Load("libAliL3Src.so"); + gSystem->Load("libAliL3Misc.so"); + gSystem->Load("libAliL3Hough.so"); + gSystem->Load("libAliL3Comp.so"); + } + } + + AliReconstructor* reconstructor = NULL; + // first check if a plugin is defined for the reconstructor + TPluginHandler* pluginHandler = + pluginManager->FindHandler("AliReconstructor", detName); + // if not, add a plugin for it + if (!pluginHandler) { + AliDebug(1, Form("defining plugin for %s", recName.Data())); + TString libs = gSystem->GetLibraries(); + if (libs.Contains("lib" + detName + "base.so") || + (gSystem->Load("lib" + detName + "base.so") >= 0)) { + pluginManager->AddHandler("AliReconstructor", detName, + recName, detName + "rec", recName + "()"); + } else { + pluginManager->AddHandler("AliReconstructor", detName, + recName, detName, recName + "()"); + } + pluginHandler = pluginManager->FindHandler("AliReconstructor", detName); + } + if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) { + reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0); + } + if (reconstructor) { + TObject* obj = fOptions.FindObject(detName.Data()); + if (obj) reconstructor->SetOption(obj->GetTitle()); + reconstructor->Init(fRunLoader); + fReconstructor[iDet] = reconstructor; + } + + // get or create the loader + if (detName != "HLT") { + fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader"); + if (!fLoader[iDet]) { + AliConfig::Instance() + ->CreateDetectorFolders(fRunLoader->GetEventFolder(), + detName, detName); + // first check if a plugin is defined for the loader + TPluginHandler* pluginHandler = + pluginManager->FindHandler("AliLoader", detName); + // if not, add a plugin for it + if (!pluginHandler) { + TString loaderName = "Ali" + detName + "Loader"; + AliDebug(1, Form("defining plugin for %s", loaderName.Data())); + pluginManager->AddHandler("AliLoader", detName, + loaderName, detName + "base", + loaderName + "(const char*, TFolder*)"); + pluginHandler = pluginManager->FindHandler("AliLoader", detName); + } + if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) { + fLoader[iDet] = + (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), + fRunLoader->GetEventFolder()); + } + if (!fLoader[iDet]) { // use default loader + fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder()); + } + if (!fLoader[iDet]) { + AliWarning(Form("couldn't get loader for %s", detName.Data())); + if (fStopOnError) return NULL; + } else { + fRunLoader->AddLoader(fLoader[iDet]); + fRunLoader->CdGAFile(); + if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE"); + fRunLoader->Write(0, TObject::kOverwrite); + } + } + } + + return reconstructor; +} + +//_____________________________________________________________________________ +Bool_t AliReconstruction::CreateVertexer() +{ +// create the vertexer + + fVertexer = NULL; + AliReconstructor* itsReconstructor = GetReconstructor(0); + if (itsReconstructor) { + fVertexer = itsReconstructor->CreateVertexer(fRunLoader); + } + if (!fVertexer) { + AliWarning("couldn't create a vertexer for ITS"); + if (fStopOnError) return kFALSE; + } + + return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliReconstruction::CreateTrackers(const TString& detectors) +{ +// create the trackers + + TString detStr = detectors; + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + if (!IsSelected(fgkDetectorName[iDet], detStr)) continue; + AliReconstructor* reconstructor = GetReconstructor(iDet); + if (!reconstructor) continue; + TString detName = fgkDetectorName[iDet]; + if (detName == "HLT") { + fRunHLTTracking = kTRUE; + continue; + } + + fTracker[iDet] = reconstructor->CreateTracker(fRunLoader); + if (!fTracker[iDet] && (iDet < 7)) { + AliWarning(Form("couldn't create a tracker for %s", detName.Data())); + if (fStopOnError) return kFALSE; + } + } + + return kTRUE; +} + +//_____________________________________________________________________________ +void AliReconstruction::CleanUp(TFile* file, TFile* fileOld) +{ +// delete trackers and the run loader and close and delete the file + + for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) { + delete fReconstructor[iDet]; + fReconstructor[iDet] = NULL; + fLoader[iDet] = NULL; + delete fTracker[iDet]; + fTracker[iDet] = NULL; + } + delete fVertexer; + fVertexer = NULL; + + delete fRunLoader; + fRunLoader = NULL; + delete fRawReader; + fRawReader = NULL; + + if (file) { + file->Close(); + delete file; + } + + if (fileOld) { + fileOld->Close(); + delete fileOld; + gSystem->Unlink("AliESDs.old.root"); + } +} + + +//_____________________________________________________________________________ +Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const +{ +// read the ESD event from a file + + if (!esd) return kFALSE; + char fileName[256]; + sprintf(fileName, "ESD_%d.%d_%s.root", + esd->GetRunNumber(), esd->GetEventNumber(), recStep); + if (gSystem->AccessPathName(fileName)) return kFALSE; + + AliDebug(1, Form("reading ESD from file %s", fileName)); + TFile* file = TFile::Open(fileName); + if (!file || !file->IsOpen()) { + AliError(Form("opening %s failed", fileName)); + delete file; + return kFALSE; + } + + gROOT->cd(); + delete esd; + esd = (AliESD*) file->Get("ESD"); + file->Close(); + delete file; + return kTRUE; +} + +//_____________________________________________________________________________ +void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const +{ +// write the ESD event to a file + + if (!esd) return; + char fileName[256]; + sprintf(fileName, "ESD_%d.%d_%s.root", + esd->GetRunNumber(), esd->GetEventNumber(), recStep); + + AliDebug(1, Form("writing ESD to file %s", fileName)); + TFile* file = TFile::Open(fileName, "recreate"); + if (!file || !file->IsOpen()) { + AliError(Form("opening %s failed", fileName)); + } else { + esd->Write("ESD"); + file->Close(); + } + delete file; +}