1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // class for running the reconstruction //
22 // Clusters and tracks are created for all detectors and all events by //
25 // AliReconstruction rec; //
28 // The Run method returns kTRUE in case of successful execution. //
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 //
33 // rec.SetInput("..."); //
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 //
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 //
45 // rec.SetEventRange(..., ...); //
47 // The index -1 (default) can be used for the last event to indicate no //
48 // upper limit of the event range. //
50 // The name of the galice file can be changed from the default //
51 // "galice.root" by passing it as argument to the AliReconstruction //
52 // constructor or by //
54 // rec.SetGAliceFile("..."); //
56 // The local reconstruction can be switched on or off for individual //
59 // rec.SetRunLocalReconstruction("..."); //
61 // The argument is a (case sensitive) string with the names of the //
62 // detectors separated by a space. The special string "ALL" selects all //
63 // available detectors. This is the default. //
65 // The reconstruction of the primary vertex position can be switched off by //
67 // rec.SetRunVertexFinder(kFALSE); //
69 // The tracking and the creation of ESD tracks can be switched on for //
70 // selected detectors by //
72 // rec.SetRunTracking("..."); //
74 // Uniform/nonuniform field tracking switches (default: uniform field) //
76 // rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
78 // The filling of additional ESD information can be steered by //
80 // rec.SetFillESD("..."); //
82 // Again, for both methods the string specifies the list of detectors. //
83 // The default is "ALL". //
85 // The call of the shortcut method //
87 // rec.SetRunReconstruction("..."); //
89 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and //
90 // SetFillESD with the same detector selecting string as argument. //
92 // The reconstruction requires digits or raw data as input. For the creation //
93 // of digits and raw data have a look at the class AliSimulation. //
95 // For debug purposes the method SetCheckPointLevel can be used. If the //
96 // argument is greater than 0, files with ESD events will be written after //
97 // selected steps of the reconstruction for each event: //
98 // level 1: after tracking and after filling of ESD (final) //
99 // level 2: in addition after each tracking step //
100 // level 3: in addition after the filling of ESD for each detector //
101 // If a final check point file exists for an event, this event will be //
102 // skipped in the reconstruction. The tracking and the filling of ESD for //
103 // a detector will be skipped as well, if the corresponding check point //
104 // file exists. The ESD event will then be loaded from the file instead. //
106 ///////////////////////////////////////////////////////////////////////////////
112 #include <TPluginManager.h>
113 #include <TStopwatch.h>
114 #include <TGeoManager.h>
115 #include <TLorentzVector.h>
117 #include "AliReconstruction.h"
118 #include "AliReconstructor.h"
120 #include "AliRunLoader.h"
122 #include "AliRawReaderFile.h"
123 #include "AliRawReaderDate.h"
124 #include "AliRawReaderRoot.h"
126 #include "AliESDfriend.h"
127 #include "AliESDVertex.h"
128 #include "AliMultiplicity.h"
129 #include "AliTracker.h"
130 #include "AliVertexer.h"
131 #include "AliVertexerTracks.h"
132 #include "AliHeader.h"
133 #include "AliGenEventHeader.h"
135 #include "AliESDpid.h"
136 #include "AliESDtrack.h"
138 #include "AliRunTag.h"
139 #include "AliDetectorTag.h"
140 #include "AliEventTag.h"
142 #include "AliTrackPointArray.h"
143 #include "AliCDBManager.h"
144 #include "AliCDBEntry.h"
145 #include "AliAlignObj.h"
147 #include "AliCentralTrigger.h"
148 #include "AliCTPRawStream.h"
150 ClassImp(AliReconstruction)
153 //_____________________________________________________________________________
154 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
156 //_____________________________________________________________________________
157 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
158 const char* name, const char* title) :
161 fUniformField(kTRUE),
162 fRunVertexFinder(kTRUE),
163 fRunHLTTracking(kFALSE),
164 fStopOnError(kFALSE),
165 fWriteAlignmentData(kFALSE),
166 fWriteESDfriend(kFALSE),
167 fFillTriggerESD(kTRUE),
169 fRunLocalReconstruction("ALL"),
172 fGAliceFileName(gAliceFilename),
179 fLoadAlignFromCDB(kTRUE),
180 fLoadAlignData("ALL"),
187 fAlignObjArray(NULL),
191 // create reconstruction object with default parameters
193 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
194 fReconstructor[iDet] = NULL;
195 fLoader[iDet] = NULL;
196 fTracker[iDet] = NULL;
201 //_____________________________________________________________________________
202 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
205 fUniformField(rec.fUniformField),
206 fRunVertexFinder(rec.fRunVertexFinder),
207 fRunHLTTracking(rec.fRunHLTTracking),
208 fStopOnError(rec.fStopOnError),
209 fWriteAlignmentData(rec.fWriteAlignmentData),
210 fWriteESDfriend(rec.fWriteESDfriend),
211 fFillTriggerESD(rec.fFillTriggerESD),
213 fRunLocalReconstruction(rec.fRunLocalReconstruction),
214 fRunTracking(rec.fRunTracking),
215 fFillESD(rec.fFillESD),
216 fGAliceFileName(rec.fGAliceFileName),
218 fEquipIdMap(rec.fEquipIdMap),
219 fFirstEvent(rec.fFirstEvent),
220 fLastEvent(rec.fLastEvent),
223 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
224 fLoadAlignData(rec.fLoadAlignData),
231 fAlignObjArray(rec.fAlignObjArray),
232 fCDBUri(rec.fCDBUri),
237 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
238 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
240 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
241 fReconstructor[iDet] = NULL;
242 fLoader[iDet] = NULL;
243 fTracker[iDet] = NULL;
245 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
246 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
250 //_____________________________________________________________________________
251 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
253 // assignment operator
255 this->~AliReconstruction();
256 new(this) AliReconstruction(rec);
260 //_____________________________________________________________________________
261 AliReconstruction::~AliReconstruction()
267 fSpecCDBUri.Delete();
270 //_____________________________________________________________________________
271 void AliReconstruction::InitCDBStorage()
273 // activate a default CDB storage
274 // First check if we have any CDB storage set, because it is used
275 // to retrieve the calibration and alignment constants
277 AliCDBManager* man = AliCDBManager::Instance();
278 if (man->IsDefaultStorageSet())
280 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
281 AliWarning("Default CDB storage has been already set !");
282 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
283 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
287 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
288 AliWarning(Form("Default CDB storage is set to: %s",fCDBUri.Data()));
289 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
290 man->SetDefaultStorage(fCDBUri);
293 // Now activate the detector specific CDB storage locations
294 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
295 TString detName = fgkDetectorName[iDet];
296 TObject* obj = fSpecCDBUri.FindObject(detName.Data());
298 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
299 AliWarning(Form("Specific CDB storage for %s is set to: %s",detName.Data(),obj->GetTitle()));
300 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
301 man->SetSpecificStorage(detName.Data(),obj->GetTitle());
307 //_____________________________________________________________________________
308 void AliReconstruction::SetDefaultStorage(const char* uri) {
309 // Store the desired default CDB storage location
310 // Activate it later within the Run() method
316 //_____________________________________________________________________________
317 void AliReconstruction::SetSpecificStorage(const char* detName, const char* uri) {
318 // Store a detector-specific CDB storage location
319 // Activate it later within the Run() method
321 TObject* obj = fSpecCDBUri.FindObject(detName);
322 if (obj) fSpecCDBUri.Remove(obj);
323 fSpecCDBUri.Add(new TNamed(detName, uri));
327 //_____________________________________________________________________________
328 Bool_t AliReconstruction::SetRunNumber()
330 // The method is called in Run() in order
331 // to set a correct run number.
332 // In case of raw data reconstruction the
333 // run number is taken from the raw data header
335 if(AliCDBManager::Instance()->GetRun() < 0) {
337 AliError("No run loader is found !");
340 // read run number from gAlice
341 if(fRunLoader->GetAliRun())
342 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
345 if(fRawReader->NextEvent()) {
346 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
347 fRawReader->RewindEvents();
350 AliError("No raw-data events found !");
355 AliError("Neither gAlice nor RawReader objects are found !");
359 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
364 //_____________________________________________________________________________
365 Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
367 // Read collection of alignment objects (AliAlignObj derived) saved
368 // in the TClonesArray ClArrayName and apply them to the geometry
369 // manager singleton.
372 Int_t nvols = alObjArray->GetEntriesFast();
376 for(Int_t j=0; j<nvols; j++)
378 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
379 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
382 if (AliDebugLevelClass() >= 1) {
383 gGeoManager->GetTopNode()->CheckOverlaps(20);
384 TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
385 if(ovexlist->GetEntriesFast()){
386 AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
394 //_____________________________________________________________________________
395 Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
397 // Fills array of single detector's alignable objects from CDB
399 AliDebug(2, Form("Loading alignment data for detector: %s",detName));
403 AliCDBPath path(detName,"Align","Data");
405 entry=AliCDBManager::Instance()->Get(path.GetPath());
407 AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
411 TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
412 alignArray->SetOwner(0);
413 AliDebug(2,Form("Found %d alignment objects for %s",
414 alignArray->GetEntries(),detName));
416 AliAlignObj *alignObj=0;
417 TIter iter(alignArray);
419 // loop over align objects in detector
420 while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
421 fAlignObjArray->Add(alignObj);
423 // delete entry --- Don't delete, it is cached!
425 AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
430 //_____________________________________________________________________________
431 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
433 // Read the alignment objects from CDB.
434 // Each detector is supposed to have the
435 // alignment objects in DET/Align/Data CDB path.
436 // All the detector objects are then collected,
437 // sorted by geometry level (starting from ALIC) and
438 // then applied to the TGeo geometry.
439 // Finally an overlaps check is performed.
441 // Load alignment data from CDB and fill fAlignObjArray
442 if(fLoadAlignFromCDB){
443 if(!fAlignObjArray) fAlignObjArray = new TObjArray();
445 //fAlignObjArray->RemoveAll();
446 fAlignObjArray->Clear();
447 fAlignObjArray->SetOwner(0);
449 TString detStr = detectors;
450 TString dataNotLoaded="";
451 TString dataLoaded="";
453 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
454 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
455 if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
456 dataNotLoaded += fgkDetectorName[iDet];
457 dataNotLoaded += " ";
459 dataLoaded += fgkDetectorName[iDet];
462 } // end loop over detectors
464 if ((detStr.CompareTo("ALL") == 0)) detStr = "";
465 dataNotLoaded += detStr;
466 AliInfo(Form("Alignment data loaded for: %s",
468 AliInfo(Form("Didn't/couldn't load alignment data for: %s",
469 dataNotLoaded.Data()));
470 } // fLoadAlignFromCDB flag
472 // Check if the array with alignment objects was
473 // provided by the user. If yes, apply the objects
474 // to the present TGeo geometry
475 if (fAlignObjArray) {
476 if (gGeoManager && gGeoManager->IsClosed()) {
477 if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
478 AliError("The misalignment of one or more volumes failed!"
479 "Compare the list of simulated detectors and the list of detector alignment data!");
484 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
492 //_____________________________________________________________________________
493 void AliReconstruction::SetGAliceFile(const char* fileName)
495 // set the name of the galice file
497 fGAliceFileName = fileName;
500 //_____________________________________________________________________________
501 void AliReconstruction::SetOption(const char* detector, const char* option)
503 // set options for the reconstruction of a detector
505 TObject* obj = fOptions.FindObject(detector);
506 if (obj) fOptions.Remove(obj);
507 fOptions.Add(new TNamed(detector, option));
511 //_____________________________________________________________________________
512 Bool_t AliReconstruction::Run(const char* input,
513 Int_t firstEvent, Int_t lastEvent)
515 // run the reconstruction
518 if (!input) input = fInput.Data();
519 TString fileName(input);
520 if (fileName.EndsWith("/")) {
521 fRawReader = new AliRawReaderFile(fileName);
522 } else if (fileName.EndsWith(".root")) {
523 fRawReader = new AliRawReaderRoot(fileName);
524 } else if (!fileName.IsNull()) {
525 fRawReader = new AliRawReaderDate(fileName);
526 fRawReader->SelectEvents(7);
528 if (!fEquipIdMap.IsNull() && fRawReader)
529 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
532 // get the run loader
533 if (!InitRunLoader()) return kFALSE;
535 // Initialize the CDB storage
538 // Set run number in CDBManager (if it is not already set by the user)
539 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
541 // Import ideal TGeo geometry and apply misalignment
543 TString geom(gSystem->DirName(fGAliceFileName));
544 geom += "/geometry.root";
545 TGeoManager::Import(geom.Data());
546 if (!gGeoManager) if (fStopOnError) return kFALSE;
548 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
550 // Temporary fix by A.Gheata
551 // Could be removed with the next Root version (>5.11)
553 TIter next(gGeoManager->GetListOfVolumes());
555 while ((vol = (TGeoVolume *)next())) {
556 if (vol->GetVoxels()) {
557 if (vol->GetVoxels()->NeedRebuild()) {
558 vol->GetVoxels()->Voxelize();
565 // local reconstruction
566 if (!fRunLocalReconstruction.IsNull()) {
567 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
568 if (fStopOnError) {CleanUp(); return kFALSE;}
571 // if (!fRunVertexFinder && fRunTracking.IsNull() &&
572 // fFillESD.IsNull()) return kTRUE;
575 if (fRunVertexFinder && !CreateVertexer()) {
583 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
591 TStopwatch stopwatch;
594 // get the possibly already existing ESD file and tree
595 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
596 TFile* fileOld = NULL;
597 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
598 if (!gSystem->AccessPathName("AliESDs.root")){
599 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
600 fileOld = TFile::Open("AliESDs.old.root");
601 if (fileOld && fileOld->IsOpen()) {
602 treeOld = (TTree*) fileOld->Get("esdTree");
603 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
604 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
605 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
609 // create the ESD output file and tree
610 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
611 if (!file->IsOpen()) {
612 AliError("opening AliESDs.root failed");
613 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
615 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
616 tree->Branch("ESD", "AliESD", &esd);
617 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
618 hlttree->Branch("ESD", "AliESD", &hltesd);
619 delete esd; delete hltesd;
620 esd = NULL; hltesd = NULL;
622 // create the file and tree with ESD additions
623 TFile *filef=0; TTree *treef=0; AliESDfriend *esdf=0;
624 if (fWriteESDfriend) {
625 filef = TFile::Open("AliESDfriends.root", "RECREATE");
626 if (!filef->IsOpen()) {
627 AliError("opening AliESDfriends.root failed");
629 treef = new TTree("esdFriendTree", "Tree with ESD friends");
630 treef->Branch("ESDfriend", "AliESDfriend", &esdf);
635 AliVertexerTracks tVertexer;
638 if (fRawReader) fRawReader->RewindEvents();
640 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
641 if (fRawReader) fRawReader->NextEvent();
642 if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
643 // copy old ESD to the new one
645 treeOld->SetBranchAddress("ESD", &esd);
646 treeOld->GetEntry(iEvent);
650 hlttreeOld->SetBranchAddress("ESD", &hltesd);
651 hlttreeOld->GetEntry(iEvent);
657 AliInfo(Form("processing event %d", iEvent));
658 fRunLoader->GetEvent(iEvent);
661 sprintf(fileName, "ESD_%d.%d_final.root",
662 fRunLoader->GetHeader()->GetRun(),
663 fRunLoader->GetHeader()->GetEventNrInRun());
664 if (!gSystem->AccessPathName(fileName)) continue;
666 // local reconstruction
667 if (!fRunLocalReconstruction.IsNull()) {
668 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
669 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
673 esd = new AliESD; hltesd = new AliESD;
674 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
675 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
676 esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
677 hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
679 // Set magnetic field from the tracker
680 esd->SetMagneticField(AliTracker::GetBz());
681 hltesd->SetMagneticField(AliTracker::GetBz());
684 if (fRunVertexFinder) {
685 if (!ReadESD(esd, "vertex")) {
686 if (!RunVertexFinder(esd)) {
687 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
689 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
694 if (!fRunTracking.IsNull()) {
695 if (fRunHLTTracking) {
696 hltesd->SetVertex(esd->GetVertex());
697 if (!RunHLTTracking(hltesd)) {
698 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
704 if (!fRunTracking.IsNull()) {
705 if (!ReadESD(esd, "tracking")) {
706 if (!RunTracking(esd)) {
707 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
709 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
714 if (!fFillESD.IsNull()) {
715 if (!FillESD(esd, fFillESD)) {
716 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
721 AliESDpid::MakePID(esd);
722 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
724 if (fFillTriggerESD) {
725 if (!ReadESD(esd, "trigger")) {
726 if (!FillTriggerESD(esd)) {
727 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
729 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
733 esd->SetPrimaryVertex(tVertexer.FindPrimaryVertex(esd));
741 if (fWriteESDfriend) {
742 esdf=new AliESDfriend();
743 esd->GetESDfriend(esdf);
747 if (fCheckPointLevel > 0) WriteESD(esd, "final");
749 delete esd; delete hltesd;
750 esd = NULL; hltesd = NULL;
753 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
754 stopwatch.RealTime(),stopwatch.CpuTime()));
760 if (fWriteESDfriend) {
762 treef->Write(); delete treef; filef->Close(); delete filef;
765 // Create tags for the events in the ESD tree (the ESD tree is always present)
766 // In case of empty events the tags will contain dummy values
768 CleanUp(file, fileOld);
774 //_____________________________________________________________________________
775 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
777 // run the local reconstruction
779 TStopwatch stopwatch;
782 TString detStr = detectors;
783 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
784 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
785 AliReconstructor* reconstructor = GetReconstructor(iDet);
786 if (!reconstructor) continue;
787 if (reconstructor->HasLocalReconstruction()) continue;
789 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
790 TStopwatch stopwatchDet;
791 stopwatchDet.Start();
793 fRawReader->RewindEvents();
794 reconstructor->Reconstruct(fRunLoader, fRawReader);
796 reconstructor->Reconstruct(fRunLoader);
798 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
799 fgkDetectorName[iDet],
800 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
803 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
804 AliError(Form("the following detectors were not found: %s",
806 if (fStopOnError) return kFALSE;
809 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
810 stopwatch.RealTime(),stopwatch.CpuTime()));
815 //_____________________________________________________________________________
816 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
818 // run the local reconstruction
820 TStopwatch stopwatch;
823 TString detStr = detectors;
824 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
825 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
826 AliReconstructor* reconstructor = GetReconstructor(iDet);
827 if (!reconstructor) continue;
828 AliLoader* loader = fLoader[iDet];
830 // conversion of digits
831 if (fRawReader && reconstructor->HasDigitConversion()) {
832 AliInfo(Form("converting raw data digits into root objects for %s",
833 fgkDetectorName[iDet]));
834 TStopwatch stopwatchDet;
835 stopwatchDet.Start();
836 loader->LoadDigits("update");
837 loader->CleanDigits();
838 loader->MakeDigitsContainer();
839 TTree* digitsTree = loader->TreeD();
840 reconstructor->ConvertDigits(fRawReader, digitsTree);
841 loader->WriteDigits("OVERWRITE");
842 loader->UnloadDigits();
843 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
844 fgkDetectorName[iDet],
845 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
848 // local reconstruction
849 if (!reconstructor->HasLocalReconstruction()) continue;
850 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
851 TStopwatch stopwatchDet;
852 stopwatchDet.Start();
853 loader->LoadRecPoints("update");
854 loader->CleanRecPoints();
855 loader->MakeRecPointsContainer();
856 TTree* clustersTree = loader->TreeR();
857 if (fRawReader && !reconstructor->HasDigitConversion()) {
858 reconstructor->Reconstruct(fRawReader, clustersTree);
860 loader->LoadDigits("read");
861 TTree* digitsTree = loader->TreeD();
863 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
864 if (fStopOnError) return kFALSE;
866 reconstructor->Reconstruct(digitsTree, clustersTree);
868 loader->UnloadDigits();
870 loader->WriteRecPoints("OVERWRITE");
871 loader->UnloadRecPoints();
872 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
873 fgkDetectorName[iDet],
874 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
877 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
878 AliError(Form("the following detectors were not found: %s",
880 if (fStopOnError) return kFALSE;
883 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
884 stopwatch.RealTime(),stopwatch.CpuTime()));
889 //_____________________________________________________________________________
890 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
892 // run the barrel tracking
894 TStopwatch stopwatch;
897 AliESDVertex* vertex = NULL;
898 Double_t vtxPos[3] = {0, 0, 0};
899 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
901 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
902 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
903 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
907 AliInfo("running the ITS vertex finder");
908 if (fLoader[0]) fLoader[0]->LoadRecPoints();
909 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
910 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
912 AliWarning("Vertex not found");
913 vertex = new AliESDVertex();
914 vertex->SetName("default");
917 vertex->SetTruePos(vtxPos); // store also the vertex from MC
918 vertex->SetName("reconstructed");
922 AliInfo("getting the primary vertex from MC");
923 vertex = new AliESDVertex(vtxPos, vtxErr);
927 vertex->GetXYZ(vtxPos);
928 vertex->GetSigmaXYZ(vtxErr);
930 AliWarning("no vertex reconstructed");
931 vertex = new AliESDVertex(vtxPos, vtxErr);
933 esd->SetVertex(vertex);
934 // if SPD multiplicity has been determined, it is stored in the ESD
935 AliMultiplicity *mult= fVertexer->GetMultiplicity();
936 if(mult)esd->SetMultiplicity(mult);
938 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
939 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
943 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
944 stopwatch.RealTime(),stopwatch.CpuTime()));
949 //_____________________________________________________________________________
950 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
952 // run the HLT barrel tracking
954 TStopwatch stopwatch;
958 AliError("Missing runLoader!");
962 AliInfo("running HLT tracking");
964 // Get a pointer to the HLT reconstructor
965 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
966 if (!reconstructor) return kFALSE;
969 for (Int_t iDet = 1; iDet >= 0; iDet--) {
970 TString detName = fgkDetectorName[iDet];
971 AliDebug(1, Form("%s HLT tracking", detName.Data()));
972 reconstructor->SetOption(detName.Data());
973 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
975 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
976 if (fStopOnError) return kFALSE;
980 Double_t vtxErr[3]={0.005,0.005,0.010};
981 const AliESDVertex *vertex = esd->GetVertex();
982 vertex->GetXYZ(vtxPos);
983 tracker->SetVertex(vtxPos,vtxErr);
985 fLoader[iDet]->LoadRecPoints("read");
986 TTree* tree = fLoader[iDet]->TreeR();
988 AliError(Form("Can't get the %s cluster tree", detName.Data()));
991 tracker->LoadClusters(tree);
993 if (tracker->Clusters2Tracks(esd) != 0) {
994 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
998 tracker->UnloadClusters();
1003 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1004 stopwatch.RealTime(),stopwatch.CpuTime()));
1009 //_____________________________________________________________________________
1010 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1012 // run the barrel tracking
1014 TStopwatch stopwatch;
1017 AliInfo("running tracking");
1019 // pass 1: TPC + ITS inwards
1020 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1021 if (!fTracker[iDet]) continue;
1022 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1025 fLoader[iDet]->LoadRecPoints("read");
1026 TTree* tree = fLoader[iDet]->TreeR();
1028 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1031 fTracker[iDet]->LoadClusters(tree);
1034 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1035 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1038 if (fCheckPointLevel > 1) {
1039 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1041 // preliminary PID in TPC needed by the ITS tracker
1043 GetReconstructor(1)->FillESD(fRunLoader, esd);
1044 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1045 AliESDpid::MakePID(esd);
1049 // pass 2: ALL backwards
1050 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1051 if (!fTracker[iDet]) continue;
1052 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1055 if (iDet > 1) { // all except ITS, TPC
1057 fLoader[iDet]->LoadRecPoints("read");
1058 tree = fLoader[iDet]->TreeR();
1060 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1063 fTracker[iDet]->LoadClusters(tree);
1067 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1068 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1071 if (fCheckPointLevel > 1) {
1072 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1076 if (iDet > 2) { // all except ITS, TPC, TRD
1077 fTracker[iDet]->UnloadClusters();
1078 fLoader[iDet]->UnloadRecPoints();
1080 // updated PID in TPC needed by the ITS tracker -MI
1082 GetReconstructor(1)->FillESD(fRunLoader, esd);
1083 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1084 AliESDpid::MakePID(esd);
1088 // write space-points to the ESD in case alignment data output
1090 if (fWriteAlignmentData)
1091 WriteAlignmentData(esd);
1093 // pass 3: TRD + TPC + ITS refit inwards
1094 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1095 if (!fTracker[iDet]) continue;
1096 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1099 if (fTracker[iDet]->RefitInward(esd) != 0) {
1100 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1103 if (fCheckPointLevel > 1) {
1104 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1108 fTracker[iDet]->UnloadClusters();
1109 fLoader[iDet]->UnloadRecPoints();
1112 // Propagate track to the vertex - if not done by ITS
1114 Int_t ntracks = esd->GetNumberOfTracks();
1115 for (Int_t itrack=0; itrack<ntracks; itrack++){
1116 const Double_t kRadius = 3; // beam pipe radius
1117 const Double_t kMaxStep = 5; // max step
1118 const Double_t kMaxD = 123456; // max distance to prim vertex
1119 Double_t fieldZ = AliTracker::GetBz(); //
1120 AliESDtrack * track = esd->GetTrack(itrack);
1121 if (!track) continue;
1122 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1123 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1124 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1127 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1128 stopwatch.RealTime(),stopwatch.CpuTime()));
1133 //_____________________________________________________________________________
1134 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1136 // fill the event summary data
1138 TStopwatch stopwatch;
1140 AliInfo("filling ESD");
1142 TString detStr = detectors;
1143 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1144 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1145 AliReconstructor* reconstructor = GetReconstructor(iDet);
1146 if (!reconstructor) continue;
1148 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1149 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1150 TTree* clustersTree = NULL;
1151 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1152 fLoader[iDet]->LoadRecPoints("read");
1153 clustersTree = fLoader[iDet]->TreeR();
1154 if (!clustersTree) {
1155 AliError(Form("Can't get the %s clusters tree",
1156 fgkDetectorName[iDet]));
1157 if (fStopOnError) return kFALSE;
1160 if (fRawReader && !reconstructor->HasDigitConversion()) {
1161 reconstructor->FillESD(fRawReader, clustersTree, esd);
1163 TTree* digitsTree = NULL;
1164 if (fLoader[iDet]) {
1165 fLoader[iDet]->LoadDigits("read");
1166 digitsTree = fLoader[iDet]->TreeD();
1168 AliError(Form("Can't get the %s digits tree",
1169 fgkDetectorName[iDet]));
1170 if (fStopOnError) return kFALSE;
1173 reconstructor->FillESD(digitsTree, clustersTree, esd);
1174 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1176 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1177 fLoader[iDet]->UnloadRecPoints();
1181 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1183 reconstructor->FillESD(fRunLoader, esd);
1185 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1189 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1190 AliError(Form("the following detectors were not found: %s",
1192 if (fStopOnError) return kFALSE;
1195 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1196 stopwatch.RealTime(),stopwatch.CpuTime()));
1201 //_____________________________________________________________________________
1202 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1204 // Reads the trigger decision which is
1205 // stored in Trigger.root file and fills
1206 // the corresponding esd entries
1208 AliInfo("Filling trigger information into the ESD");
1211 AliCTPRawStream input(fRawReader);
1212 if (!input.Next()) {
1213 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1216 esd->SetTriggerMask(input.GetClassMask());
1217 esd->SetTriggerCluster(input.GetClusterMask());
1220 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1222 if (!runloader->LoadTrigger()) {
1223 AliCentralTrigger *aCTP = runloader->GetTrigger();
1224 esd->SetTriggerMask(aCTP->GetClassMask());
1225 esd->SetTriggerCluster(aCTP->GetClusterMask());
1228 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1233 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1241 //_____________________________________________________________________________
1242 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1244 // check whether detName is contained in detectors
1245 // if yes, it is removed from detectors
1247 // check if all detectors are selected
1248 if ((detectors.CompareTo("ALL") == 0) ||
1249 detectors.BeginsWith("ALL ") ||
1250 detectors.EndsWith(" ALL") ||
1251 detectors.Contains(" ALL ")) {
1256 // search for the given detector
1257 Bool_t result = kFALSE;
1258 if ((detectors.CompareTo(detName) == 0) ||
1259 detectors.BeginsWith(detName+" ") ||
1260 detectors.EndsWith(" "+detName) ||
1261 detectors.Contains(" "+detName+" ")) {
1262 detectors.ReplaceAll(detName, "");
1266 // clean up the detectors string
1267 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1268 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1269 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1274 //_____________________________________________________________________________
1275 Bool_t AliReconstruction::InitRunLoader()
1277 // get or create the run loader
1279 if (gAlice) delete gAlice;
1282 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1283 // load all base libraries to get the loader classes
1284 TString libs = gSystem->GetLibraries();
1285 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1286 TString detName = fgkDetectorName[iDet];
1287 if (detName == "HLT") continue;
1288 if (libs.Contains("lib" + detName + "base.so")) continue;
1289 gSystem->Load("lib" + detName + "base.so");
1291 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1293 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1297 fRunLoader->CdGAFile();
1298 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1299 if (fRunLoader->LoadgAlice() == 0) {
1300 gAlice = fRunLoader->GetAliRun();
1301 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1304 if (!gAlice && !fRawReader) {
1305 AliError(Form("no gAlice object found in file %s",
1306 fGAliceFileName.Data()));
1311 } else { // galice.root does not exist
1313 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1317 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1318 AliConfig::GetDefaultEventFolderName(),
1321 AliError(Form("could not create run loader in file %s",
1322 fGAliceFileName.Data()));
1326 fRunLoader->MakeTree("E");
1328 while (fRawReader->NextEvent()) {
1329 fRunLoader->SetEventNumber(iEvent);
1330 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1332 fRunLoader->MakeTree("H");
1333 fRunLoader->TreeE()->Fill();
1336 fRawReader->RewindEvents();
1337 fRunLoader->WriteHeader("OVERWRITE");
1338 fRunLoader->CdGAFile();
1339 fRunLoader->Write(0, TObject::kOverwrite);
1340 // AliTracker::SetFieldMap(???);
1346 //_____________________________________________________________________________
1347 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1349 // get the reconstructor object and the loader for a detector
1351 if (fReconstructor[iDet]) return fReconstructor[iDet];
1353 // load the reconstructor object
1354 TPluginManager* pluginManager = gROOT->GetPluginManager();
1355 TString detName = fgkDetectorName[iDet];
1356 TString recName = "Ali" + detName + "Reconstructor";
1357 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1359 if (detName == "HLT") {
1360 if (!gROOT->GetClass("AliLevel3")) {
1361 gSystem->Load("libAliL3Src.so");
1362 gSystem->Load("libAliL3Misc.so");
1363 gSystem->Load("libAliL3Hough.so");
1364 gSystem->Load("libAliL3Comp.so");
1368 AliReconstructor* reconstructor = NULL;
1369 // first check if a plugin is defined for the reconstructor
1370 TPluginHandler* pluginHandler =
1371 pluginManager->FindHandler("AliReconstructor", detName);
1372 // if not, add a plugin for it
1373 if (!pluginHandler) {
1374 AliDebug(1, Form("defining plugin for %s", recName.Data()));
1375 TString libs = gSystem->GetLibraries();
1376 if (libs.Contains("lib" + detName + "base.so") ||
1377 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1378 pluginManager->AddHandler("AliReconstructor", detName,
1379 recName, detName + "rec", recName + "()");
1381 pluginManager->AddHandler("AliReconstructor", detName,
1382 recName, detName, recName + "()");
1384 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1386 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1387 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1389 if (reconstructor) {
1390 TObject* obj = fOptions.FindObject(detName.Data());
1391 if (obj) reconstructor->SetOption(obj->GetTitle());
1392 reconstructor->Init(fRunLoader);
1393 fReconstructor[iDet] = reconstructor;
1396 // get or create the loader
1397 if (detName != "HLT") {
1398 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1399 if (!fLoader[iDet]) {
1400 AliConfig::Instance()
1401 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1403 // first check if a plugin is defined for the loader
1404 TPluginHandler* pluginHandler =
1405 pluginManager->FindHandler("AliLoader", detName);
1406 // if not, add a plugin for it
1407 if (!pluginHandler) {
1408 TString loaderName = "Ali" + detName + "Loader";
1409 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1410 pluginManager->AddHandler("AliLoader", detName,
1411 loaderName, detName + "base",
1412 loaderName + "(const char*, TFolder*)");
1413 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1415 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1417 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1418 fRunLoader->GetEventFolder());
1420 if (!fLoader[iDet]) { // use default loader
1421 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1423 if (!fLoader[iDet]) {
1424 AliWarning(Form("couldn't get loader for %s", detName.Data()));
1425 if (fStopOnError) return NULL;
1427 fRunLoader->AddLoader(fLoader[iDet]);
1428 fRunLoader->CdGAFile();
1429 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1430 fRunLoader->Write(0, TObject::kOverwrite);
1435 return reconstructor;
1438 //_____________________________________________________________________________
1439 Bool_t AliReconstruction::CreateVertexer()
1441 // create the vertexer
1444 AliReconstructor* itsReconstructor = GetReconstructor(0);
1445 if (itsReconstructor) {
1446 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1449 AliWarning("couldn't create a vertexer for ITS");
1450 if (fStopOnError) return kFALSE;
1456 //_____________________________________________________________________________
1457 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1459 // create the trackers
1461 TString detStr = detectors;
1462 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1463 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1464 AliReconstructor* reconstructor = GetReconstructor(iDet);
1465 if (!reconstructor) continue;
1466 TString detName = fgkDetectorName[iDet];
1467 if (detName == "HLT") {
1468 fRunHLTTracking = kTRUE;
1472 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1473 if (!fTracker[iDet] && (iDet < 7)) {
1474 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1475 if (fStopOnError) return kFALSE;
1482 //_____________________________________________________________________________
1483 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1485 // delete trackers and the run loader and close and delete the file
1487 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1488 delete fReconstructor[iDet];
1489 fReconstructor[iDet] = NULL;
1490 fLoader[iDet] = NULL;
1491 delete fTracker[iDet];
1492 fTracker[iDet] = NULL;
1510 gSystem->Unlink("AliESDs.old.root");
1515 //_____________________________________________________________________________
1516 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1518 // read the ESD event from a file
1520 if (!esd) return kFALSE;
1522 sprintf(fileName, "ESD_%d.%d_%s.root",
1523 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1524 if (gSystem->AccessPathName(fileName)) return kFALSE;
1526 AliInfo(Form("reading ESD from file %s", fileName));
1527 AliDebug(1, Form("reading ESD from file %s", fileName));
1528 TFile* file = TFile::Open(fileName);
1529 if (!file || !file->IsOpen()) {
1530 AliError(Form("opening %s failed", fileName));
1537 esd = (AliESD*) file->Get("ESD");
1543 //_____________________________________________________________________________
1544 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1546 // write the ESD event to a file
1550 sprintf(fileName, "ESD_%d.%d_%s.root",
1551 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1553 AliDebug(1, Form("writing ESD to file %s", fileName));
1554 TFile* file = TFile::Open(fileName, "recreate");
1555 if (!file || !file->IsOpen()) {
1556 AliError(Form("opening %s failed", fileName));
1567 //_____________________________________________________________________________
1568 void AliReconstruction::CreateTag(TFile* file)
1573 Double_t fMUONMASS = 0.105658369;
1576 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1577 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1579 TLorentzVector fEPvector;
1581 Float_t fZVertexCut = 10.0;
1582 Float_t fRhoVertexCut = 2.0;
1584 Float_t fLowPtCut = 1.0;
1585 Float_t fHighPtCut = 3.0;
1586 Float_t fVeryHighPtCut = 10.0;
1589 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1591 // Creates the tags for all the events in a given ESD file
1593 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1594 Int_t nPos, nNeg, nNeutr;
1595 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1596 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1597 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1598 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1599 Float_t maxPt = .0, meanPt = .0, totalP = .0;
1601 Int_t iRunNumber = 0;
1602 TString fVertexName("default");
1604 AliRunTag *tag = new AliRunTag();
1605 AliEventTag *evTag = new AliEventTag();
1606 TTree ttag("T","A Tree with event tags");
1607 TBranch * btag = ttag.Branch("AliTAG", &tag);
1608 btag->SetCompressionLevel(9);
1610 AliInfo(Form("Creating the tags......."));
1612 if (!file || !file->IsOpen()) {
1613 AliError(Form("opening failed"));
1617 Int_t firstEvent = 0,lastEvent = 0;
1618 TTree *t = (TTree*) file->Get("esdTree");
1619 TBranch * b = t->GetBranch("ESD");
1621 b->SetAddress(&esd);
1624 Int_t iInitRunNumber = esd->GetRunNumber();
1626 Int_t iNumberOfEvents = b->GetEntries();
1627 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1655 b->GetEntry(iEventNumber);
1656 iRunNumber = esd->GetRunNumber();
1657 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1658 const AliESDVertex * vertexIn = esd->GetVertex();
1659 if (!vertexIn) AliError("ESD has not defined vertex.");
1660 if (vertexIn) fVertexName = vertexIn->GetName();
1661 if(fVertexName != "default") fVertexflag = 1;
1662 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1663 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1664 UInt_t status = esdTrack->GetStatus();
1666 //select only tracks with ITS refit
1667 if ((status&AliESDtrack::kITSrefit)==0) continue;
1668 //select only tracks with TPC refit
1669 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1671 //select only tracks with the "combined PID"
1672 if ((status&AliESDtrack::kESDpid)==0) continue;
1674 esdTrack->GetPxPyPz(p);
1675 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1676 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1679 if(fPt > maxPt) maxPt = fPt;
1681 if(esdTrack->GetSign() > 0) {
1683 if(fPt > fLowPtCut) nCh1GeV++;
1684 if(fPt > fHighPtCut) nCh3GeV++;
1685 if(fPt > fVeryHighPtCut) nCh10GeV++;
1687 if(esdTrack->GetSign() < 0) {
1689 if(fPt > fLowPtCut) nCh1GeV++;
1690 if(fPt > fHighPtCut) nCh3GeV++;
1691 if(fPt > fVeryHighPtCut) nCh10GeV++;
1693 if(esdTrack->GetSign() == 0) nNeutr++;
1697 esdTrack->GetESDpid(prob);
1700 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1701 if(rcc == 0.0) continue;
1704 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1707 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1709 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1711 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1713 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1715 if(fPt > fLowPtCut) nEl1GeV++;
1716 if(fPt > fHighPtCut) nEl3GeV++;
1717 if(fPt > fVeryHighPtCut) nEl10GeV++;
1725 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1726 // loop over all reconstructed tracks (also first track of combination)
1727 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1728 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1729 if (muonTrack == 0x0) continue;
1731 // Coordinates at vertex
1732 fZ = muonTrack->GetZ();
1733 fY = muonTrack->GetBendingCoor();
1734 fX = muonTrack->GetNonBendingCoor();
1736 fThetaX = muonTrack->GetThetaX();
1737 fThetaY = muonTrack->GetThetaY();
1739 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1740 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1741 fPxRec = fPzRec * TMath::Tan(fThetaX);
1742 fPyRec = fPzRec * TMath::Tan(fThetaY);
1743 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1745 //ChiSquare of the track if needed
1746 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1747 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1748 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1750 // total number of muons inside a vertex cut
1751 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1753 if(fEPvector.Pt() > fLowPtCut) {
1755 if(fEPvector.Pt() > fHighPtCut) {
1757 if (fEPvector.Pt() > fVeryHighPtCut) {
1765 // Fill the event tags
1767 meanPt = meanPt/ntrack;
1769 evTag->SetEventId(iEventNumber+1);
1771 evTag->SetVertexX(vertexIn->GetXv());
1772 evTag->SetVertexY(vertexIn->GetYv());
1773 evTag->SetVertexZ(vertexIn->GetZv());
1774 evTag->SetVertexZError(vertexIn->GetZRes());
1776 evTag->SetVertexFlag(fVertexflag);
1778 evTag->SetT0VertexZ(esd->GetT0zVertex());
1780 evTag->SetTriggerMask(esd->GetTriggerMask());
1781 evTag->SetTriggerCluster(esd->GetTriggerCluster());
1783 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1784 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1785 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1786 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1787 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1788 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1791 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1792 evTag->SetNumOfPosTracks(nPos);
1793 evTag->SetNumOfNegTracks(nNeg);
1794 evTag->SetNumOfNeutrTracks(nNeutr);
1796 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1797 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1798 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1799 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1801 evTag->SetNumOfProtons(nProtons);
1802 evTag->SetNumOfKaons(nKaons);
1803 evTag->SetNumOfPions(nPions);
1804 evTag->SetNumOfMuons(nMuons);
1805 evTag->SetNumOfElectrons(nElectrons);
1806 evTag->SetNumOfPhotons(nGamas);
1807 evTag->SetNumOfPi0s(nPi0s);
1808 evTag->SetNumOfNeutrons(nNeutrons);
1809 evTag->SetNumOfKaon0s(nK0s);
1811 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1812 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1813 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1814 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1815 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1816 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1817 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1818 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1819 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1821 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1822 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1824 evTag->SetTotalMomentum(totalP);
1825 evTag->SetMeanPt(meanPt);
1826 evTag->SetMaxPt(maxPt);
1828 tag->SetRunId(iInitRunNumber);
1829 tag->AddEventTag(*evTag);
1831 lastEvent = iNumberOfEvents;
1837 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1838 tag->GetRunId(),firstEvent,lastEvent );
1839 AliInfo(Form("writing tags to file %s", fileName));
1840 AliDebug(1, Form("writing tags to file %s", fileName));
1842 TFile* ftag = TFile::Open(fileName, "recreate");
1851 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1853 // Write space-points which are then used in the alignment procedures
1854 // For the moment only ITS, TRD and TPC
1856 // Load TOF clusters
1858 fLoader[3]->LoadRecPoints("read");
1859 TTree* tree = fLoader[3]->TreeR();
1861 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1864 fTracker[3]->LoadClusters(tree);
1866 Int_t ntracks = esd->GetNumberOfTracks();
1867 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1869 AliESDtrack *track = esd->GetTrack(itrack);
1872 for (Int_t iDet = 3; iDet >= 0; iDet--)
1873 nsp += track->GetNcls(iDet);
1875 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1876 track->SetTrackPointArray(sp);
1878 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1879 AliTracker *tracker = fTracker[iDet];
1880 if (!tracker) continue;
1881 Int_t nspdet = track->GetNcls(iDet);
1882 if (nspdet <= 0) continue;
1883 track->GetClusters(iDet,idx);
1887 while (isp < nspdet) {
1888 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1889 const Int_t kNTPCmax = 159;
1890 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
1891 if (!isvalid) continue;
1892 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1898 fTracker[3]->UnloadClusters();
1899 fLoader[3]->UnloadRecPoints();