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 // In case of raw-data reconstruction the user can modify the default //
51 // number of events per digits/clusters/tracks file. In case the option //
52 // is not used the number is set 1. In case the user provides 0, than //
53 // the number of events is equal to the number of events inside the //
54 // raw-data file (i.e. one digits/clusters/tracks file): //
56 // rec.SetNumberOfEventsPerFile(...); //
59 // The name of the galice file can be changed from the default //
60 // "galice.root" by passing it as argument to the AliReconstruction //
61 // constructor or by //
63 // rec.SetGAliceFile("..."); //
65 // The local reconstruction can be switched on or off for individual //
68 // rec.SetRunLocalReconstruction("..."); //
70 // The argument is a (case sensitive) string with the names of the //
71 // detectors separated by a space. The special string "ALL" selects all //
72 // available detectors. This is the default. //
74 // The reconstruction of the primary vertex position can be switched off by //
76 // rec.SetRunVertexFinder(kFALSE); //
78 // The tracking and the creation of ESD tracks can be switched on for //
79 // selected detectors by //
81 // rec.SetRunTracking("..."); //
83 // Uniform/nonuniform field tracking switches (default: uniform field) //
85 // rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
87 // The filling of additional ESD information can be steered by //
89 // rec.SetFillESD("..."); //
91 // Again, for both methods the string specifies the list of detectors. //
92 // The default is "ALL". //
94 // The call of the shortcut method //
96 // rec.SetRunReconstruction("..."); //
98 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and //
99 // SetFillESD with the same detector selecting string as argument. //
101 // The reconstruction requires digits or raw data as input. For the creation //
102 // of digits and raw data have a look at the class AliSimulation. //
104 // For debug purposes the method SetCheckPointLevel can be used. If the //
105 // argument is greater than 0, files with ESD events will be written after //
106 // selected steps of the reconstruction for each event: //
107 // level 1: after tracking and after filling of ESD (final) //
108 // level 2: in addition after each tracking step //
109 // level 3: in addition after the filling of ESD for each detector //
110 // If a final check point file exists for an event, this event will be //
111 // skipped in the reconstruction. The tracking and the filling of ESD for //
112 // a detector will be skipped as well, if the corresponding check point //
113 // file exists. The ESD event will then be loaded from the file instead. //
115 ///////////////////////////////////////////////////////////////////////////////
121 #include <TPluginManager.h>
122 #include <TStopwatch.h>
123 #include <TGeoManager.h>
124 #include <TLorentzVector.h>
126 #include "AliReconstruction.h"
127 #include "AliReconstructor.h"
129 #include "AliRunLoader.h"
131 #include "AliRawReaderFile.h"
132 #include "AliRawReaderDate.h"
133 #include "AliRawReaderRoot.h"
134 #include "AliRawEventHeaderBase.h"
136 #include "AliESDfriend.h"
137 #include "AliESDVertex.h"
138 #include "AliMultiplicity.h"
139 #include "AliTracker.h"
140 #include "AliVertexer.h"
141 #include "AliVertexerTracks.h"
142 #include "AliV0vertexer.h"
143 #include "AliCascadeVertexer.h"
144 #include "AliHeader.h"
145 #include "AliGenEventHeader.h"
147 #include "AliESDpid.h"
148 #include "AliESDtrack.h"
150 #include "AliRunTag.h"
151 #include "AliDetectorTag.h"
152 #include "AliEventTag.h"
154 #include "AliTrackPointArray.h"
155 #include "AliCDBManager.h"
156 #include "AliCDBEntry.h"
157 #include "AliAlignObj.h"
159 #include "AliCentralTrigger.h"
160 #include "AliCTPRawStream.h"
162 ClassImp(AliReconstruction)
165 //_____________________________________________________________________________
166 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
168 //_____________________________________________________________________________
169 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
170 const char* name, const char* title) :
173 fUniformField(kTRUE),
174 fRunVertexFinder(kTRUE),
175 fRunHLTTracking(kFALSE),
176 fRunMuonTracking(kFALSE),
177 fStopOnError(kFALSE),
178 fWriteAlignmentData(kFALSE),
179 fWriteESDfriend(kFALSE),
181 fFillTriggerESD(kTRUE),
183 fRunLocalReconstruction("ALL"),
186 fGAliceFileName(gAliceFilename),
191 fNumberOfEventsPerFile(1),
194 fLoadAlignFromCDB(kTRUE),
195 fLoadAlignData("ALL"),
201 fDiamondProfile(NULL),
203 fAlignObjArray(NULL),
207 // create reconstruction object with default parameters
209 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
210 fReconstructor[iDet] = NULL;
211 fLoader[iDet] = NULL;
212 fTracker[iDet] = NULL;
217 //_____________________________________________________________________________
218 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
221 fUniformField(rec.fUniformField),
222 fRunVertexFinder(rec.fRunVertexFinder),
223 fRunHLTTracking(rec.fRunHLTTracking),
224 fRunMuonTracking(rec.fRunMuonTracking),
225 fStopOnError(rec.fStopOnError),
226 fWriteAlignmentData(rec.fWriteAlignmentData),
227 fWriteESDfriend(rec.fWriteESDfriend),
228 fWriteAOD(rec.fWriteAOD),
229 fFillTriggerESD(rec.fFillTriggerESD),
231 fRunLocalReconstruction(rec.fRunLocalReconstruction),
232 fRunTracking(rec.fRunTracking),
233 fFillESD(rec.fFillESD),
234 fGAliceFileName(rec.fGAliceFileName),
236 fEquipIdMap(rec.fEquipIdMap),
237 fFirstEvent(rec.fFirstEvent),
238 fLastEvent(rec.fLastEvent),
239 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
242 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
243 fLoadAlignData(rec.fLoadAlignData),
249 fDiamondProfile(NULL),
251 fAlignObjArray(rec.fAlignObjArray),
252 fCDBUri(rec.fCDBUri),
257 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
258 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
260 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
261 fReconstructor[iDet] = NULL;
262 fLoader[iDet] = NULL;
263 fTracker[iDet] = NULL;
265 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
266 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
270 //_____________________________________________________________________________
271 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
273 // assignment operator
275 this->~AliReconstruction();
276 new(this) AliReconstruction(rec);
280 //_____________________________________________________________________________
281 AliReconstruction::~AliReconstruction()
287 fSpecCDBUri.Delete();
290 //_____________________________________________________________________________
291 void AliReconstruction::InitCDBStorage()
293 // activate a default CDB storage
294 // First check if we have any CDB storage set, because it is used
295 // to retrieve the calibration and alignment constants
297 AliCDBManager* man = AliCDBManager::Instance();
298 if (man->IsDefaultStorageSet())
300 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
301 AliWarning("Default CDB storage has been already set !");
302 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
303 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
307 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
308 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
309 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
310 man->SetDefaultStorage(fCDBUri);
313 // Now activate the detector specific CDB storage locations
314 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
315 TObject* obj = fSpecCDBUri[i];
317 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
318 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
319 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
320 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
325 //_____________________________________________________________________________
326 void AliReconstruction::SetDefaultStorage(const char* uri) {
327 // Store the desired default CDB storage location
328 // Activate it later within the Run() method
334 //_____________________________________________________________________________
335 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
336 // Store a detector-specific CDB storage location
337 // Activate it later within the Run() method
339 AliCDBPath aPath(calibType);
340 if(!aPath.IsValid()){
341 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
342 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
343 if(!strcmp(calibType, fgkDetectorName[iDet])) {
344 aPath.SetPath(Form("%s/*", calibType));
345 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
349 if(!aPath.IsValid()){
350 AliError(Form("Not a valid path or detector: %s", calibType));
355 // check that calibType refers to a "valid" detector name
356 Bool_t isDetector = kFALSE;
357 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
358 TString detName = fgkDetectorName[iDet];
359 if(aPath.GetLevel0() == detName) {
366 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
370 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
371 if (obj) fSpecCDBUri.Remove(obj);
372 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
376 //_____________________________________________________________________________
377 Bool_t AliReconstruction::SetRunNumber()
379 // The method is called in Run() in order
380 // to set a correct run number.
381 // In case of raw data reconstruction the
382 // run number is taken from the raw data header
384 if(AliCDBManager::Instance()->GetRun() < 0) {
386 AliError("No run loader is found !");
389 // read run number from gAlice
390 if(fRunLoader->GetAliRun())
391 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
394 if(fRawReader->NextEvent()) {
395 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
396 fRawReader->RewindEvents();
399 AliError("No raw-data events found !");
404 AliError("Neither gAlice nor RawReader objects are found !");
408 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
413 //_____________________________________________________________________________
414 Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
416 // Read collection of alignment objects (AliAlignObj derived) saved
417 // in the TClonesArray ClArrayName and apply them to the geometry
418 // manager singleton.
421 Int_t nvols = alObjArray->GetEntriesFast();
425 for(Int_t j=0; j<nvols; j++)
427 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
428 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
431 if (AliDebugLevelClass() >= 1) {
432 gGeoManager->GetTopNode()->CheckOverlaps(20);
433 TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
434 if(ovexlist->GetEntriesFast()){
435 AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
443 //_____________________________________________________________________________
444 Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
446 // Fills array of single detector's alignable objects from CDB
448 AliDebug(2, Form("Loading alignment data for detector: %s",detName));
452 AliCDBPath path(detName,"Align","Data");
454 entry=AliCDBManager::Instance()->Get(path.GetPath());
456 AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
460 TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
461 alignArray->SetOwner(0);
462 AliDebug(2,Form("Found %d alignment objects for %s",
463 alignArray->GetEntries(),detName));
465 AliAlignObj *alignObj=0;
466 TIter iter(alignArray);
468 // loop over align objects in detector
469 while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
470 fAlignObjArray->Add(alignObj);
472 // delete entry --- Don't delete, it is cached!
474 AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
479 //_____________________________________________________________________________
480 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
482 // Read the alignment objects from CDB.
483 // Each detector is supposed to have the
484 // alignment objects in DET/Align/Data CDB path.
485 // All the detector objects are then collected,
486 // sorted by geometry level (starting from ALIC) and
487 // then applied to the TGeo geometry.
488 // Finally an overlaps check is performed.
490 // Load alignment data from CDB and fill fAlignObjArray
491 if(fLoadAlignFromCDB){
492 if(!fAlignObjArray) fAlignObjArray = new TObjArray();
494 //fAlignObjArray->RemoveAll();
495 fAlignObjArray->Clear();
496 fAlignObjArray->SetOwner(0);
498 TString detStr = detectors;
499 TString dataNotLoaded="";
500 TString dataLoaded="";
502 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
503 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
504 if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
505 dataNotLoaded += fgkDetectorName[iDet];
506 dataNotLoaded += " ";
508 dataLoaded += fgkDetectorName[iDet];
511 } // end loop over detectors
513 if ((detStr.CompareTo("ALL") == 0)) detStr = "";
514 dataNotLoaded += detStr;
515 if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
517 if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
518 dataNotLoaded.Data()));
519 } // fLoadAlignFromCDB flag
521 // Check if the array with alignment objects was
522 // provided by the user. If yes, apply the objects
523 // to the present TGeo geometry
524 if (fAlignObjArray) {
525 if (gGeoManager && gGeoManager->IsClosed()) {
526 if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
527 AliError("The misalignment of one or more volumes failed!"
528 "Compare the list of simulated detectors and the list of detector alignment data!");
533 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
538 delete fAlignObjArray; fAlignObjArray=0;
540 // Update the TGeoPhysicalNodes
541 gGeoManager->RefreshPhysicalNodes();
546 //_____________________________________________________________________________
547 void AliReconstruction::SetGAliceFile(const char* fileName)
549 // set the name of the galice file
551 fGAliceFileName = fileName;
554 //_____________________________________________________________________________
555 void AliReconstruction::SetOption(const char* detector, const char* option)
557 // set options for the reconstruction of a detector
559 TObject* obj = fOptions.FindObject(detector);
560 if (obj) fOptions.Remove(obj);
561 fOptions.Add(new TNamed(detector, option));
565 //_____________________________________________________________________________
566 Bool_t AliReconstruction::Run(const char* input)
568 // run the reconstruction
571 if (!input) input = fInput.Data();
572 TString fileName(input);
573 if (fileName.EndsWith("/")) {
574 fRawReader = new AliRawReaderFile(fileName);
575 } else if (fileName.EndsWith(".root")) {
576 fRawReader = new AliRawReaderRoot(fileName);
577 } else if (!fileName.IsNull()) {
578 fRawReader = new AliRawReaderDate(fileName);
579 fRawReader->SelectEvents(7);
581 if (!fEquipIdMap.IsNull() && fRawReader)
582 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
585 // get the run loader
586 if (!InitRunLoader()) return kFALSE;
588 // Initialize the CDB storage
591 // Set run number in CDBManager (if it is not already set by the user)
592 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
594 // Import ideal TGeo geometry and apply misalignment
596 TString geom(gSystem->DirName(fGAliceFileName));
597 geom += "/geometry.root";
598 TGeoManager::Import(geom.Data());
599 if (!gGeoManager) if (fStopOnError) return kFALSE;
602 AliCDBManager* man = AliCDBManager::Instance();
603 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
605 // local reconstruction
606 if (!fRunLocalReconstruction.IsNull()) {
607 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
608 if (fStopOnError) {CleanUp(); return kFALSE;}
611 // if (!fRunVertexFinder && fRunTracking.IsNull() &&
612 // fFillESD.IsNull()) return kTRUE;
615 if (fRunVertexFinder && !CreateVertexer()) {
623 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
631 TStopwatch stopwatch;
634 // get the possibly already existing ESD file and tree
635 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
636 TFile* fileOld = NULL;
637 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
638 if (!gSystem->AccessPathName("AliESDs.root")){
639 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
640 fileOld = TFile::Open("AliESDs.old.root");
641 if (fileOld && fileOld->IsOpen()) {
642 treeOld = (TTree*) fileOld->Get("esdTree");
643 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
644 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
645 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
649 // create the ESD output file and tree
650 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
651 if (!file->IsOpen()) {
652 AliError("opening AliESDs.root failed");
653 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
655 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
656 tree->Branch("ESD", "AliESD", &esd);
657 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
658 hlttree->Branch("ESD", "AliESD", &hltesd);
659 delete esd; delete hltesd;
660 esd = NULL; hltesd = NULL;
662 // create the branch with ESD additions
663 AliESDfriend *esdf=0;
664 if (fWriteESDfriend) {
665 TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
666 br->SetFile("AliESDfriends.root");
669 AliVertexerTracks tVertexer;
670 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
673 if (fRawReader) fRawReader->RewindEvents();
675 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
676 if (fRawReader) fRawReader->NextEvent();
677 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
678 // copy old ESD to the new one
680 treeOld->SetBranchAddress("ESD", &esd);
681 treeOld->GetEntry(iEvent);
685 hlttreeOld->SetBranchAddress("ESD", &hltesd);
686 hlttreeOld->GetEntry(iEvent);
692 AliInfo(Form("processing event %d", iEvent));
693 fRunLoader->GetEvent(iEvent);
696 sprintf(aFileName, "ESD_%d.%d_final.root",
697 fRunLoader->GetHeader()->GetRun(),
698 fRunLoader->GetHeader()->GetEventNrInRun());
699 if (!gSystem->AccessPathName(aFileName)) continue;
701 // local reconstruction
702 if (!fRunLocalReconstruction.IsNull()) {
703 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
704 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
708 esd = new AliESD; hltesd = new AliESD;
709 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
710 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
711 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
712 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
714 // Set magnetic field from the tracker
715 esd->SetMagneticField(AliTracker::GetBz());
716 hltesd->SetMagneticField(AliTracker::GetBz());
718 // Fill raw-data error log into the ESD
719 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
722 if (fRunVertexFinder) {
723 if (!ReadESD(esd, "vertex")) {
724 if (!RunVertexFinder(esd)) {
725 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
727 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
732 if (!fRunTracking.IsNull()) {
733 if (fRunHLTTracking) {
734 hltesd->SetVertex(esd->GetVertex());
735 if (!RunHLTTracking(hltesd)) {
736 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
742 if (!fRunTracking.IsNull()) {
743 if (fRunMuonTracking) {
744 if (!RunMuonTracking()) {
745 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
751 if (!fRunTracking.IsNull()) {
752 if (!ReadESD(esd, "tracking")) {
753 if (!RunTracking(esd)) {
754 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
756 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
761 if (!fFillESD.IsNull()) {
762 if (!FillESD(esd, fFillESD)) {
763 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
766 // fill Event header information from the RawEventHeader
767 if (fRawReader){FillRawEventHeaderESD(esd);}
770 AliESDpid::MakePID(esd);
771 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
773 if (fFillTriggerESD) {
774 if (!ReadESD(esd, "trigger")) {
775 if (!FillTriggerESD(esd)) {
776 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
778 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
782 //Try to improve the reconstructed primary vertex position using the tracks
783 AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
785 if (pvtx->GetStatus()) {
786 // Store the improved primary vertex
787 esd->SetPrimaryVertex(pvtx);
788 // Propagate the tracks to the DCA to the improved primary vertex
789 Double_t somethingbig = 777.;
790 Double_t bz = esd->GetMagneticField();
791 Int_t nt=esd->GetNumberOfTracks();
793 AliESDtrack *t = esd->GetTrack(nt);
794 t->RelateToVertex(pvtx, bz, somethingbig);
801 vtxer.Tracks2V0vertices(esd);
804 AliCascadeVertexer cvtxer;
805 cvtxer.V0sTracks2CascadeVertices(esd);
809 if (fWriteESDfriend) {
810 esdf=new AliESDfriend();
811 esd->GetESDfriend(esdf);
818 if (fCheckPointLevel > 0) WriteESD(esd, "final");
820 delete esd; delete esdf; delete hltesd;
821 esd = NULL; esdf=NULL; hltesd = NULL;
824 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
825 stopwatch.RealTime(),stopwatch.CpuTime()));
829 tree->SetBranchStatus("ESDfriend*",0);
837 // Create tags for the events in the ESD tree (the ESD tree is always present)
838 // In case of empty events the tags will contain dummy values
840 CleanUp(file, fileOld);
846 //_____________________________________________________________________________
847 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
849 // run the local reconstruction
851 TStopwatch stopwatch;
854 AliCDBManager* man = AliCDBManager::Instance();
855 Bool_t origCache = man->GetCacheFlag();
857 TString detStr = detectors;
858 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
859 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
860 AliReconstructor* reconstructor = GetReconstructor(iDet);
861 if (!reconstructor) continue;
862 if (reconstructor->HasLocalReconstruction()) continue;
864 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
865 TStopwatch stopwatchDet;
866 stopwatchDet.Start();
868 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
870 man->SetCacheFlag(kTRUE);
871 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
872 man->GetAll(calibPath); // entries are cached!
875 fRawReader->RewindEvents();
876 reconstructor->Reconstruct(fRunLoader, fRawReader);
878 reconstructor->Reconstruct(fRunLoader);
880 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
881 fgkDetectorName[iDet],
882 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
884 // unload calibration data
888 man->SetCacheFlag(origCache);
890 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
891 AliError(Form("the following detectors were not found: %s",
893 if (fStopOnError) return kFALSE;
896 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
897 stopwatch.RealTime(),stopwatch.CpuTime()));
902 //_____________________________________________________________________________
903 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
905 // run the local reconstruction
907 TStopwatch stopwatch;
910 TString detStr = detectors;
911 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
912 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
913 AliReconstructor* reconstructor = GetReconstructor(iDet);
914 if (!reconstructor) continue;
915 AliLoader* loader = fLoader[iDet];
917 // conversion of digits
918 if (fRawReader && reconstructor->HasDigitConversion()) {
919 AliInfo(Form("converting raw data digits into root objects for %s",
920 fgkDetectorName[iDet]));
921 TStopwatch stopwatchDet;
922 stopwatchDet.Start();
923 loader->LoadDigits("update");
924 loader->CleanDigits();
925 loader->MakeDigitsContainer();
926 TTree* digitsTree = loader->TreeD();
927 reconstructor->ConvertDigits(fRawReader, digitsTree);
928 loader->WriteDigits("OVERWRITE");
929 loader->UnloadDigits();
930 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
931 fgkDetectorName[iDet],
932 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
935 // local reconstruction
936 if (!reconstructor->HasLocalReconstruction()) continue;
937 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
938 TStopwatch stopwatchDet;
939 stopwatchDet.Start();
940 loader->LoadRecPoints("update");
941 loader->CleanRecPoints();
942 loader->MakeRecPointsContainer();
943 TTree* clustersTree = loader->TreeR();
944 if (fRawReader && !reconstructor->HasDigitConversion()) {
945 reconstructor->Reconstruct(fRawReader, clustersTree);
947 loader->LoadDigits("read");
948 TTree* digitsTree = loader->TreeD();
950 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
951 if (fStopOnError) return kFALSE;
953 reconstructor->Reconstruct(digitsTree, clustersTree);
955 loader->UnloadDigits();
957 loader->WriteRecPoints("OVERWRITE");
958 loader->UnloadRecPoints();
959 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
960 fgkDetectorName[iDet],
961 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
964 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
965 AliError(Form("the following detectors were not found: %s",
967 if (fStopOnError) return kFALSE;
970 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
971 stopwatch.RealTime(),stopwatch.CpuTime()));
976 //_____________________________________________________________________________
977 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
979 // run the barrel tracking
981 TStopwatch stopwatch;
984 AliESDVertex* vertex = NULL;
985 Double_t vtxPos[3] = {0, 0, 0};
986 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
988 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
989 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
990 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
994 AliInfo("running the ITS vertex finder");
995 if (fLoader[0]) fLoader[0]->LoadRecPoints();
996 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
997 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
999 AliWarning("Vertex not found");
1000 vertex = new AliESDVertex();
1001 vertex->SetName("default");
1004 vertex->SetTruePos(vtxPos); // store also the vertex from MC
1005 vertex->SetName("reconstructed");
1009 AliInfo("getting the primary vertex from MC");
1010 vertex = new AliESDVertex(vtxPos, vtxErr);
1014 vertex->GetXYZ(vtxPos);
1015 vertex->GetSigmaXYZ(vtxErr);
1017 AliWarning("no vertex reconstructed");
1018 vertex = new AliESDVertex(vtxPos, vtxErr);
1020 esd->SetVertex(vertex);
1021 // if SPD multiplicity has been determined, it is stored in the ESD
1022 AliMultiplicity *mult= fVertexer->GetMultiplicity();
1023 if(mult)esd->SetMultiplicity(mult);
1025 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1026 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1030 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1031 stopwatch.RealTime(),stopwatch.CpuTime()));
1036 //_____________________________________________________________________________
1037 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1039 // run the HLT barrel tracking
1041 TStopwatch stopwatch;
1045 AliError("Missing runLoader!");
1049 AliInfo("running HLT tracking");
1051 // Get a pointer to the HLT reconstructor
1052 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1053 if (!reconstructor) return kFALSE;
1056 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1057 TString detName = fgkDetectorName[iDet];
1058 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1059 reconstructor->SetOption(detName.Data());
1060 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1062 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1063 if (fStopOnError) return kFALSE;
1067 Double_t vtxErr[3]={0.005,0.005,0.010};
1068 const AliESDVertex *vertex = esd->GetVertex();
1069 vertex->GetXYZ(vtxPos);
1070 tracker->SetVertex(vtxPos,vtxErr);
1072 fLoader[iDet]->LoadRecPoints("read");
1073 TTree* tree = fLoader[iDet]->TreeR();
1075 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1078 tracker->LoadClusters(tree);
1080 if (tracker->Clusters2Tracks(esd) != 0) {
1081 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1085 tracker->UnloadClusters();
1090 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1091 stopwatch.RealTime(),stopwatch.CpuTime()));
1096 //_____________________________________________________________________________
1097 Bool_t AliReconstruction::RunMuonTracking()
1099 // run the muon spectrometer tracking
1101 TStopwatch stopwatch;
1105 AliError("Missing runLoader!");
1108 Int_t iDet = 7; // for MUON
1110 AliInfo("is running...");
1112 // Get a pointer to the MUON reconstructor
1113 AliReconstructor *reconstructor = GetReconstructor(iDet);
1114 if (!reconstructor) return kFALSE;
1117 TString detName = fgkDetectorName[iDet];
1118 AliDebug(1, Form("%s tracking", detName.Data()));
1119 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1121 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1126 fLoader[iDet]->LoadTracks("update");
1127 fLoader[iDet]->CleanTracks();
1128 fLoader[iDet]->MakeTracksContainer();
1131 fLoader[iDet]->LoadRecPoints("read");
1133 if (!tracker->Clusters2Tracks(0x0)) {
1134 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1137 fLoader[iDet]->UnloadRecPoints();
1139 fLoader[iDet]->WriteTracks("OVERWRITE");
1140 fLoader[iDet]->UnloadTracks();
1145 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1146 stopwatch.RealTime(),stopwatch.CpuTime()));
1152 //_____________________________________________________________________________
1153 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1155 // run the barrel tracking
1157 TStopwatch stopwatch;
1160 AliInfo("running tracking");
1162 // pass 1: TPC + ITS inwards
1163 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1164 if (!fTracker[iDet]) continue;
1165 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1168 fLoader[iDet]->LoadRecPoints("read");
1169 TTree* tree = fLoader[iDet]->TreeR();
1171 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1174 fTracker[iDet]->LoadClusters(tree);
1177 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1178 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1181 if (fCheckPointLevel > 1) {
1182 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1184 // preliminary PID in TPC needed by the ITS tracker
1186 GetReconstructor(1)->FillESD(fRunLoader, esd);
1187 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1188 AliESDpid::MakePID(esd);
1192 // pass 2: ALL backwards
1193 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1194 if (!fTracker[iDet]) continue;
1195 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1198 if (iDet > 1) { // all except ITS, TPC
1200 fLoader[iDet]->LoadRecPoints("read");
1201 tree = fLoader[iDet]->TreeR();
1203 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1206 fTracker[iDet]->LoadClusters(tree);
1210 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1211 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1214 if (fCheckPointLevel > 1) {
1215 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1219 if (iDet > 2) { // all except ITS, TPC, TRD
1220 fTracker[iDet]->UnloadClusters();
1221 fLoader[iDet]->UnloadRecPoints();
1223 // updated PID in TPC needed by the ITS tracker -MI
1225 GetReconstructor(1)->FillESD(fRunLoader, esd);
1226 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1227 AliESDpid::MakePID(esd);
1231 // write space-points to the ESD in case alignment data output
1233 if (fWriteAlignmentData)
1234 WriteAlignmentData(esd);
1236 // pass 3: TRD + TPC + ITS refit inwards
1237 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1238 if (!fTracker[iDet]) continue;
1239 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1242 if (fTracker[iDet]->RefitInward(esd) != 0) {
1243 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1246 if (fCheckPointLevel > 1) {
1247 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1251 fTracker[iDet]->UnloadClusters();
1252 fLoader[iDet]->UnloadRecPoints();
1255 // Propagate track to the vertex - if not done by ITS
1257 Int_t ntracks = esd->GetNumberOfTracks();
1258 for (Int_t itrack=0; itrack<ntracks; itrack++){
1259 const Double_t kRadius = 3; // beam pipe radius
1260 const Double_t kMaxStep = 5; // max step
1261 const Double_t kMaxD = 123456; // max distance to prim vertex
1262 Double_t fieldZ = AliTracker::GetBz(); //
1263 AliESDtrack * track = esd->GetTrack(itrack);
1264 if (!track) continue;
1265 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1266 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1267 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1270 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1271 stopwatch.RealTime(),stopwatch.CpuTime()));
1276 //_____________________________________________________________________________
1277 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1279 // fill the event summary data
1281 TStopwatch stopwatch;
1283 AliInfo("filling ESD");
1285 TString detStr = detectors;
1286 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1287 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1288 AliReconstructor* reconstructor = GetReconstructor(iDet);
1289 if (!reconstructor) continue;
1291 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1292 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1293 TTree* clustersTree = NULL;
1294 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1295 fLoader[iDet]->LoadRecPoints("read");
1296 clustersTree = fLoader[iDet]->TreeR();
1297 if (!clustersTree) {
1298 AliError(Form("Can't get the %s clusters tree",
1299 fgkDetectorName[iDet]));
1300 if (fStopOnError) return kFALSE;
1303 if (fRawReader && !reconstructor->HasDigitConversion()) {
1304 reconstructor->FillESD(fRawReader, clustersTree, esd);
1306 TTree* digitsTree = NULL;
1307 if (fLoader[iDet]) {
1308 fLoader[iDet]->LoadDigits("read");
1309 digitsTree = fLoader[iDet]->TreeD();
1311 AliError(Form("Can't get the %s digits tree",
1312 fgkDetectorName[iDet]));
1313 if (fStopOnError) return kFALSE;
1316 reconstructor->FillESD(digitsTree, clustersTree, esd);
1317 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1319 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1320 fLoader[iDet]->UnloadRecPoints();
1324 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1326 reconstructor->FillESD(fRunLoader, esd);
1328 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1332 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1333 AliError(Form("the following detectors were not found: %s",
1335 if (fStopOnError) return kFALSE;
1338 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1339 stopwatch.RealTime(),stopwatch.CpuTime()));
1344 //_____________________________________________________________________________
1345 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1347 // Reads the trigger decision which is
1348 // stored in Trigger.root file and fills
1349 // the corresponding esd entries
1351 AliInfo("Filling trigger information into the ESD");
1354 AliCTPRawStream input(fRawReader);
1355 if (!input.Next()) {
1356 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1359 esd->SetTriggerMask(input.GetClassMask());
1360 esd->SetTriggerCluster(input.GetClusterMask());
1363 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1365 if (!runloader->LoadTrigger()) {
1366 AliCentralTrigger *aCTP = runloader->GetTrigger();
1367 esd->SetTriggerMask(aCTP->GetClassMask());
1368 esd->SetTriggerCluster(aCTP->GetClusterMask());
1371 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1376 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1388 //_____________________________________________________________________________
1389 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1392 // Filling information from RawReader Header
1395 AliInfo("Filling information from RawReader Header");
1396 esd->SetBunchCrossNumber(0);
1397 esd->SetOrbitNumber(0);
1398 esd->SetTimeStamp(0);
1399 esd->SetEventType(0);
1400 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1402 esd->SetBunchCrossNumber((eventHeader->GetP("Id")[0]));
1403 esd->SetOrbitNumber((eventHeader->GetP("Id")[1]));
1404 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1405 esd->SetEventType((eventHeader->Get("Type")));
1412 //_____________________________________________________________________________
1413 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1415 // check whether detName is contained in detectors
1416 // if yes, it is removed from detectors
1418 // check if all detectors are selected
1419 if ((detectors.CompareTo("ALL") == 0) ||
1420 detectors.BeginsWith("ALL ") ||
1421 detectors.EndsWith(" ALL") ||
1422 detectors.Contains(" ALL ")) {
1427 // search for the given detector
1428 Bool_t result = kFALSE;
1429 if ((detectors.CompareTo(detName) == 0) ||
1430 detectors.BeginsWith(detName+" ") ||
1431 detectors.EndsWith(" "+detName) ||
1432 detectors.Contains(" "+detName+" ")) {
1433 detectors.ReplaceAll(detName, "");
1437 // clean up the detectors string
1438 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1439 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1440 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1445 //_____________________________________________________________________________
1446 Bool_t AliReconstruction::InitRunLoader()
1448 // get or create the run loader
1450 if (gAlice) delete gAlice;
1453 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1454 // load all base libraries to get the loader classes
1455 TString libs = gSystem->GetLibraries();
1456 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1457 TString detName = fgkDetectorName[iDet];
1458 if (detName == "HLT") continue;
1459 if (libs.Contains("lib" + detName + "base.so")) continue;
1460 gSystem->Load("lib" + detName + "base.so");
1462 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1464 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1468 fRunLoader->CdGAFile();
1469 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1470 if (fRunLoader->LoadgAlice() == 0) {
1471 gAlice = fRunLoader->GetAliRun();
1472 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1475 if (!gAlice && !fRawReader) {
1476 AliError(Form("no gAlice object found in file %s",
1477 fGAliceFileName.Data()));
1482 } else { // galice.root does not exist
1484 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1488 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1489 AliConfig::GetDefaultEventFolderName(),
1492 AliError(Form("could not create run loader in file %s",
1493 fGAliceFileName.Data()));
1497 fRunLoader->MakeTree("E");
1499 while (fRawReader->NextEvent()) {
1500 fRunLoader->SetEventNumber(iEvent);
1501 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1503 fRunLoader->MakeTree("H");
1504 fRunLoader->TreeE()->Fill();
1507 fRawReader->RewindEvents();
1508 if (fNumberOfEventsPerFile > 0)
1509 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1511 fRunLoader->SetNumberOfEventsPerFile(iEvent);
1512 fRunLoader->WriteHeader("OVERWRITE");
1513 fRunLoader->CdGAFile();
1514 fRunLoader->Write(0, TObject::kOverwrite);
1515 // AliTracker::SetFieldMap(???);
1521 //_____________________________________________________________________________
1522 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1524 // get the reconstructor object and the loader for a detector
1526 if (fReconstructor[iDet]) return fReconstructor[iDet];
1528 // load the reconstructor object
1529 TPluginManager* pluginManager = gROOT->GetPluginManager();
1530 TString detName = fgkDetectorName[iDet];
1531 TString recName = "Ali" + detName + "Reconstructor";
1532 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1534 if (detName == "HLT") {
1535 if (!gROOT->GetClass("AliLevel3")) {
1536 gSystem->Load("libAliHLTSrc.so");
1537 gSystem->Load("libAliHLTMisc.so");
1538 gSystem->Load("libAliHLTHough.so");
1539 gSystem->Load("libAliHLTComp.so");
1543 AliReconstructor* reconstructor = NULL;
1544 // first check if a plugin is defined for the reconstructor
1545 TPluginHandler* pluginHandler =
1546 pluginManager->FindHandler("AliReconstructor", detName);
1547 // if not, add a plugin for it
1548 if (!pluginHandler) {
1549 AliDebug(1, Form("defining plugin for %s", recName.Data()));
1550 TString libs = gSystem->GetLibraries();
1551 if (libs.Contains("lib" + detName + "base.so") ||
1552 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1553 pluginManager->AddHandler("AliReconstructor", detName,
1554 recName, detName + "rec", recName + "()");
1556 pluginManager->AddHandler("AliReconstructor", detName,
1557 recName, detName, recName + "()");
1559 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1561 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1562 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1564 if (reconstructor) {
1565 TObject* obj = fOptions.FindObject(detName.Data());
1566 if (obj) reconstructor->SetOption(obj->GetTitle());
1567 reconstructor->Init(fRunLoader);
1568 fReconstructor[iDet] = reconstructor;
1571 // get or create the loader
1572 if (detName != "HLT") {
1573 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1574 if (!fLoader[iDet]) {
1575 AliConfig::Instance()
1576 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1578 // first check if a plugin is defined for the loader
1580 pluginManager->FindHandler("AliLoader", detName);
1581 // if not, add a plugin for it
1582 if (!pluginHandler) {
1583 TString loaderName = "Ali" + detName + "Loader";
1584 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1585 pluginManager->AddHandler("AliLoader", detName,
1586 loaderName, detName + "base",
1587 loaderName + "(const char*, TFolder*)");
1588 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1590 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1592 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1593 fRunLoader->GetEventFolder());
1595 if (!fLoader[iDet]) { // use default loader
1596 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1598 if (!fLoader[iDet]) {
1599 AliWarning(Form("couldn't get loader for %s", detName.Data()));
1600 if (fStopOnError) return NULL;
1602 fRunLoader->AddLoader(fLoader[iDet]);
1603 fRunLoader->CdGAFile();
1604 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1605 fRunLoader->Write(0, TObject::kOverwrite);
1610 return reconstructor;
1613 //_____________________________________________________________________________
1614 Bool_t AliReconstruction::CreateVertexer()
1616 // create the vertexer
1619 AliReconstructor* itsReconstructor = GetReconstructor(0);
1620 if (itsReconstructor) {
1621 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1624 AliWarning("couldn't create a vertexer for ITS");
1625 if (fStopOnError) return kFALSE;
1631 //_____________________________________________________________________________
1632 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1634 // create the trackers
1636 TString detStr = detectors;
1637 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1638 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1639 AliReconstructor* reconstructor = GetReconstructor(iDet);
1640 if (!reconstructor) continue;
1641 TString detName = fgkDetectorName[iDet];
1642 if (detName == "HLT") {
1643 fRunHLTTracking = kTRUE;
1646 if (detName == "MUON") {
1647 fRunMuonTracking = kTRUE;
1652 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1653 if (!fTracker[iDet] && (iDet < 7)) {
1654 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1655 if (fStopOnError) return kFALSE;
1662 //_____________________________________________________________________________
1663 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1665 // delete trackers and the run loader and close and delete the file
1667 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1668 delete fReconstructor[iDet];
1669 fReconstructor[iDet] = NULL;
1670 fLoader[iDet] = NULL;
1671 delete fTracker[iDet];
1672 fTracker[iDet] = NULL;
1676 delete fDiamondProfile;
1677 fDiamondProfile = NULL;
1692 gSystem->Unlink("AliESDs.old.root");
1697 //_____________________________________________________________________________
1698 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1700 // read the ESD event from a file
1702 if (!esd) return kFALSE;
1704 sprintf(fileName, "ESD_%d.%d_%s.root",
1705 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1706 if (gSystem->AccessPathName(fileName)) return kFALSE;
1708 AliInfo(Form("reading ESD from file %s", fileName));
1709 AliDebug(1, Form("reading ESD from file %s", fileName));
1710 TFile* file = TFile::Open(fileName);
1711 if (!file || !file->IsOpen()) {
1712 AliError(Form("opening %s failed", fileName));
1719 esd = (AliESD*) file->Get("ESD");
1725 //_____________________________________________________________________________
1726 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1728 // write the ESD event to a file
1732 sprintf(fileName, "ESD_%d.%d_%s.root",
1733 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1735 AliDebug(1, Form("writing ESD to file %s", fileName));
1736 TFile* file = TFile::Open(fileName, "recreate");
1737 if (!file || !file->IsOpen()) {
1738 AliError(Form("opening %s failed", fileName));
1749 //_____________________________________________________________________________
1750 void AliReconstruction::CreateTag(TFile* file)
1755 Double_t fMUONMASS = 0.105658369;
1758 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1759 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1761 TLorentzVector fEPvector;
1763 Float_t fZVertexCut = 10.0;
1764 Float_t fRhoVertexCut = 2.0;
1766 Float_t fLowPtCut = 1.0;
1767 Float_t fHighPtCut = 3.0;
1768 Float_t fVeryHighPtCut = 10.0;
1771 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1773 // Creates the tags for all the events in a given ESD file
1775 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1776 Int_t nPos, nNeg, nNeutr;
1777 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1778 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1779 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1780 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1781 Float_t maxPt = .0, meanPt = .0, totalP = .0;
1783 Int_t iRunNumber = 0;
1784 TString fVertexName("default");
1786 AliRunTag *tag = new AliRunTag();
1787 AliEventTag *evTag = new AliEventTag();
1788 TTree ttag("T","A Tree with event tags");
1789 TBranch * btag = ttag.Branch("AliTAG", &tag);
1790 btag->SetCompressionLevel(9);
1792 AliInfo(Form("Creating the tags......."));
1794 if (!file || !file->IsOpen()) {
1795 AliError(Form("opening failed"));
1799 Int_t lastEvent = 0;
1800 TTree *t = (TTree*) file->Get("esdTree");
1801 TBranch * b = t->GetBranch("ESD");
1803 b->SetAddress(&esd);
1805 b->GetEntry(fFirstEvent);
1806 Int_t iInitRunNumber = esd->GetRunNumber();
1808 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1809 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1810 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1838 b->GetEntry(iEventNumber);
1839 iRunNumber = esd->GetRunNumber();
1840 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1841 const AliESDVertex * vertexIn = esd->GetVertex();
1842 if (!vertexIn) AliError("ESD has not defined vertex.");
1843 if (vertexIn) fVertexName = vertexIn->GetName();
1844 if(fVertexName != "default") fVertexflag = 1;
1845 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1846 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1847 UInt_t status = esdTrack->GetStatus();
1849 //select only tracks with ITS refit
1850 if ((status&AliESDtrack::kITSrefit)==0) continue;
1851 //select only tracks with TPC refit
1852 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1854 //select only tracks with the "combined PID"
1855 if ((status&AliESDtrack::kESDpid)==0) continue;
1857 esdTrack->GetPxPyPz(p);
1858 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1859 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1862 if(fPt > maxPt) maxPt = fPt;
1864 if(esdTrack->GetSign() > 0) {
1866 if(fPt > fLowPtCut) nCh1GeV++;
1867 if(fPt > fHighPtCut) nCh3GeV++;
1868 if(fPt > fVeryHighPtCut) nCh10GeV++;
1870 if(esdTrack->GetSign() < 0) {
1872 if(fPt > fLowPtCut) nCh1GeV++;
1873 if(fPt > fHighPtCut) nCh3GeV++;
1874 if(fPt > fVeryHighPtCut) nCh10GeV++;
1876 if(esdTrack->GetSign() == 0) nNeutr++;
1880 esdTrack->GetESDpid(prob);
1883 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1884 if(rcc == 0.0) continue;
1887 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1890 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1892 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1894 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1896 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1898 if(fPt > fLowPtCut) nEl1GeV++;
1899 if(fPt > fHighPtCut) nEl3GeV++;
1900 if(fPt > fVeryHighPtCut) nEl10GeV++;
1908 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1909 // loop over all reconstructed tracks (also first track of combination)
1910 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1911 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1912 if (muonTrack == 0x0) continue;
1914 // Coordinates at vertex
1915 fZ = muonTrack->GetZ();
1916 fY = muonTrack->GetBendingCoor();
1917 fX = muonTrack->GetNonBendingCoor();
1919 fThetaX = muonTrack->GetThetaX();
1920 fThetaY = muonTrack->GetThetaY();
1922 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1923 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1924 fPxRec = fPzRec * TMath::Tan(fThetaX);
1925 fPyRec = fPzRec * TMath::Tan(fThetaY);
1926 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1928 //ChiSquare of the track if needed
1929 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1930 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1931 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1933 // total number of muons inside a vertex cut
1934 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1936 if(fEPvector.Pt() > fLowPtCut) {
1938 if(fEPvector.Pt() > fHighPtCut) {
1940 if (fEPvector.Pt() > fVeryHighPtCut) {
1948 // Fill the event tags
1950 meanPt = meanPt/ntrack;
1952 evTag->SetEventId(iEventNumber+1);
1954 evTag->SetVertexX(vertexIn->GetXv());
1955 evTag->SetVertexY(vertexIn->GetYv());
1956 evTag->SetVertexZ(vertexIn->GetZv());
1957 evTag->SetVertexZError(vertexIn->GetZRes());
1959 evTag->SetVertexFlag(fVertexflag);
1961 evTag->SetT0VertexZ(esd->GetT0zVertex());
1963 evTag->SetTriggerMask(esd->GetTriggerMask());
1964 evTag->SetTriggerCluster(esd->GetTriggerCluster());
1966 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1967 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1968 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1969 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1970 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1971 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1974 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1975 evTag->SetNumOfPosTracks(nPos);
1976 evTag->SetNumOfNegTracks(nNeg);
1977 evTag->SetNumOfNeutrTracks(nNeutr);
1979 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1980 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1981 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1982 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1984 evTag->SetNumOfProtons(nProtons);
1985 evTag->SetNumOfKaons(nKaons);
1986 evTag->SetNumOfPions(nPions);
1987 evTag->SetNumOfMuons(nMuons);
1988 evTag->SetNumOfElectrons(nElectrons);
1989 evTag->SetNumOfPhotons(nGamas);
1990 evTag->SetNumOfPi0s(nPi0s);
1991 evTag->SetNumOfNeutrons(nNeutrons);
1992 evTag->SetNumOfKaon0s(nK0s);
1994 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1995 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1996 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1997 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1998 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1999 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
2000 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
2001 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
2002 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
2004 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
2005 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2007 evTag->SetTotalMomentum(totalP);
2008 evTag->SetMeanPt(meanPt);
2009 evTag->SetMaxPt(maxPt);
2011 tag->SetRunId(iInitRunNumber);
2012 tag->AddEventTag(*evTag);
2014 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2015 else lastEvent = fLastEvent;
2021 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2022 tag->GetRunId(),fFirstEvent,lastEvent );
2023 AliInfo(Form("writing tags to file %s", fileName));
2024 AliDebug(1, Form("writing tags to file %s", fileName));
2026 TFile* ftag = TFile::Open(fileName, "recreate");
2035 //_____________________________________________________________________________
2036 void AliReconstruction::CreateAOD(TFile* esdFile)
2038 // do nothing for now
2044 void AliReconstruction::WriteAlignmentData(AliESD* esd)
2046 // Write space-points which are then used in the alignment procedures
2047 // For the moment only ITS, TRD and TPC
2049 // Load TOF clusters
2051 fLoader[3]->LoadRecPoints("read");
2052 TTree* tree = fLoader[3]->TreeR();
2054 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2057 fTracker[3]->LoadClusters(tree);
2059 Int_t ntracks = esd->GetNumberOfTracks();
2060 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2062 AliESDtrack *track = esd->GetTrack(itrack);
2065 for (Int_t iDet = 3; iDet >= 0; iDet--)
2066 nsp += track->GetNcls(iDet);
2068 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2069 track->SetTrackPointArray(sp);
2071 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2072 AliTracker *tracker = fTracker[iDet];
2073 if (!tracker) continue;
2074 Int_t nspdet = track->GetNcls(iDet);
2075 if (nspdet <= 0) continue;
2076 track->GetClusters(iDet,idx);
2080 while (isp < nspdet) {
2081 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2082 const Int_t kNTPCmax = 159;
2083 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
2084 if (!isvalid) continue;
2085 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2091 fTracker[3]->UnloadClusters();
2092 fLoader[3]->UnloadRecPoints();
2096 //_____________________________________________________________________________
2097 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2099 // The method reads the raw-data error log
2100 // accumulated within the rawReader.
2101 // It extracts the raw-data errors related to
2102 // the current event and stores them into
2103 // a TClonesArray inside the esd object.
2105 if (!fRawReader) return;
2107 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2109 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2111 if (iEvent != log->GetEventNumber()) continue;
2113 esd->AddRawDataErrorLog(log);