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 "AliGeomManager.h"
155 #include "AliTrackPointArray.h"
156 #include "AliCDBManager.h"
157 #include "AliCDBEntry.h"
158 #include "AliAlignObj.h"
160 #include "AliCentralTrigger.h"
161 #include "AliCTPRawStream.h"
163 #include "AliAODEvent.h"
164 #include "AliAODHeader.h"
165 #include "AliAODTrack.h"
166 #include "AliAODVertex.h"
167 #include "AliAODCluster.h"
170 ClassImp(AliReconstruction)
173 //_____________________________________________________________________________
174 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
176 //_____________________________________________________________________________
177 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
178 const char* name, const char* title) :
181 fUniformField(kTRUE),
182 fRunVertexFinder(kTRUE),
183 fRunHLTTracking(kFALSE),
184 fRunMuonTracking(kFALSE),
185 fStopOnError(kFALSE),
186 fWriteAlignmentData(kFALSE),
187 fWriteESDfriend(kFALSE),
189 fFillTriggerESD(kTRUE),
191 fRunLocalReconstruction("ALL"),
194 fGAliceFileName(gAliceFilename),
199 fNumberOfEventsPerFile(1),
202 fLoadAlignFromCDB(kTRUE),
203 fLoadAlignData("ALL"),
210 fDiamondProfile(NULL),
212 fAlignObjArray(NULL),
216 // create reconstruction object with default parameters
218 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
219 fReconstructor[iDet] = NULL;
220 fLoader[iDet] = NULL;
221 fTracker[iDet] = NULL;
226 //_____________________________________________________________________________
227 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
230 fUniformField(rec.fUniformField),
231 fRunVertexFinder(rec.fRunVertexFinder),
232 fRunHLTTracking(rec.fRunHLTTracking),
233 fRunMuonTracking(rec.fRunMuonTracking),
234 fStopOnError(rec.fStopOnError),
235 fWriteAlignmentData(rec.fWriteAlignmentData),
236 fWriteESDfriend(rec.fWriteESDfriend),
237 fWriteAOD(rec.fWriteAOD),
238 fFillTriggerESD(rec.fFillTriggerESD),
240 fRunLocalReconstruction(rec.fRunLocalReconstruction),
241 fRunTracking(rec.fRunTracking),
242 fFillESD(rec.fFillESD),
243 fGAliceFileName(rec.fGAliceFileName),
245 fEquipIdMap(rec.fEquipIdMap),
246 fFirstEvent(rec.fFirstEvent),
247 fLastEvent(rec.fLastEvent),
248 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
251 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
252 fLoadAlignData(rec.fLoadAlignData),
253 fESDPar(rec.fESDPar),
259 fDiamondProfile(NULL),
261 fAlignObjArray(rec.fAlignObjArray),
262 fCDBUri(rec.fCDBUri),
267 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
268 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
270 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
271 fReconstructor[iDet] = NULL;
272 fLoader[iDet] = NULL;
273 fTracker[iDet] = NULL;
275 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
276 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
280 //_____________________________________________________________________________
281 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
283 // assignment operator
285 this->~AliReconstruction();
286 new(this) AliReconstruction(rec);
290 //_____________________________________________________________________________
291 AliReconstruction::~AliReconstruction()
297 fSpecCDBUri.Delete();
300 //_____________________________________________________________________________
301 void AliReconstruction::InitCDBStorage()
303 // activate a default CDB storage
304 // First check if we have any CDB storage set, because it is used
305 // to retrieve the calibration and alignment constants
307 AliCDBManager* man = AliCDBManager::Instance();
308 if (man->IsDefaultStorageSet())
310 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
311 AliWarning("Default CDB storage has been already set !");
312 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
313 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
317 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
318 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
319 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
320 man->SetDefaultStorage(fCDBUri);
323 // Now activate the detector specific CDB storage locations
324 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
325 TObject* obj = fSpecCDBUri[i];
327 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
328 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
329 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
330 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
335 //_____________________________________________________________________________
336 void AliReconstruction::SetDefaultStorage(const char* uri) {
337 // Store the desired default CDB storage location
338 // Activate it later within the Run() method
344 //_____________________________________________________________________________
345 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
346 // Store a detector-specific CDB storage location
347 // Activate it later within the Run() method
349 AliCDBPath aPath(calibType);
350 if(!aPath.IsValid()){
351 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
352 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
353 if(!strcmp(calibType, fgkDetectorName[iDet])) {
354 aPath.SetPath(Form("%s/*", calibType));
355 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
359 if(!aPath.IsValid()){
360 AliError(Form("Not a valid path or detector: %s", calibType));
365 // check that calibType refers to a "valid" detector name
366 Bool_t isDetector = kFALSE;
367 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
368 TString detName = fgkDetectorName[iDet];
369 if(aPath.GetLevel0() == detName) {
376 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
380 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
381 if (obj) fSpecCDBUri.Remove(obj);
382 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
389 //_____________________________________________________________________________
390 Bool_t AliReconstruction::SetRunNumber()
392 // The method is called in Run() in order
393 // to set a correct run number.
394 // In case of raw data reconstruction the
395 // run number is taken from the raw data header
397 if(AliCDBManager::Instance()->GetRun() < 0) {
399 AliError("No run loader is found !");
402 // read run number from gAlice
403 if(fRunLoader->GetAliRun())
404 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
407 if(fRawReader->NextEvent()) {
408 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
409 fRawReader->RewindEvents();
412 AliError("No raw-data events found !");
417 AliError("Neither gAlice nor RawReader objects are found !");
421 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
426 //_____________________________________________________________________________
427 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
429 // Read the alignment objects from CDB.
430 // Each detector is supposed to have the
431 // alignment objects in DET/Align/Data CDB path.
432 // All the detector objects are then collected,
433 // sorted by geometry level (starting from ALIC) and
434 // then applied to the TGeo geometry.
435 // Finally an overlaps check is performed.
437 // Load alignment data from CDB and fill fAlignObjArray
438 if(fLoadAlignFromCDB){
440 TString detStr = detectors;
441 TString loadAlObjsListOfDets = "";
443 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
444 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
445 loadAlObjsListOfDets += fgkDetectorName[iDet];
446 loadAlObjsListOfDets += " ";
447 } // end loop over detectors
448 (AliGeomManager::Instance())->ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
450 // Check if the array with alignment objects was
451 // provided by the user. If yes, apply the objects
452 // to the present TGeo geometry
453 if (fAlignObjArray) {
454 if (gGeoManager && gGeoManager->IsClosed()) {
455 if ((AliGeomManager::Instance())->ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
456 AliError("The misalignment of one or more volumes failed!"
457 "Compare the list of simulated detectors and the list of detector alignment data!");
462 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
468 delete fAlignObjArray; fAlignObjArray=0;
470 // Update the TGeoPhysicalNodes
471 gGeoManager->RefreshPhysicalNodes();
476 //_____________________________________________________________________________
477 void AliReconstruction::SetGAliceFile(const char* fileName)
479 // set the name of the galice file
481 fGAliceFileName = fileName;
484 //_____________________________________________________________________________
485 void AliReconstruction::SetOption(const char* detector, const char* option)
487 // set options for the reconstruction of a detector
489 TObject* obj = fOptions.FindObject(detector);
490 if (obj) fOptions.Remove(obj);
491 fOptions.Add(new TNamed(detector, option));
495 //_____________________________________________________________________________
496 Bool_t AliReconstruction::Run(const char* input)
498 // run the reconstruction
501 if (!input) input = fInput.Data();
502 TString fileName(input);
503 if (fileName.EndsWith("/")) {
504 fRawReader = new AliRawReaderFile(fileName);
505 } else if (fileName.EndsWith(".root")) {
506 fRawReader = new AliRawReaderRoot(fileName);
507 } else if (!fileName.IsNull()) {
508 fRawReader = new AliRawReaderDate(fileName);
509 fRawReader->SelectEvents(7);
511 if (!fEquipIdMap.IsNull() && fRawReader)
512 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
515 // get the run loader
516 if (!InitRunLoader()) return kFALSE;
518 // Initialize the CDB storage
521 // Set run number in CDBManager (if it is not already set by the user)
522 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
524 // Import ideal TGeo geometry and apply misalignment
526 TString geom(gSystem->DirName(fGAliceFileName));
527 geom += "/geometry.root";
528 TGeoManager::Import(geom.Data());
529 if (!gGeoManager) if (fStopOnError) return kFALSE;
532 AliCDBManager* man = AliCDBManager::Instance();
533 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
535 // local reconstruction
536 if (!fRunLocalReconstruction.IsNull()) {
537 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
538 if (fStopOnError) {CleanUp(); return kFALSE;}
541 // if (!fRunVertexFinder && fRunTracking.IsNull() &&
542 // fFillESD.IsNull()) return kTRUE;
545 if (fRunVertexFinder && !CreateVertexer()) {
553 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
561 TStopwatch stopwatch;
564 // get the possibly already existing ESD file and tree
565 AliESD* esd = new AliESD(); AliESD* hltesd = new AliESD();
566 TFile* fileOld = NULL;
567 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
568 if (!gSystem->AccessPathName("AliESDs.root")){
569 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
570 fileOld = TFile::Open("AliESDs.old.root");
571 if (fileOld && fileOld->IsOpen()) {
572 treeOld = (TTree*) fileOld->Get("esdTree");
573 if (treeOld)esd->ReadFromTree(treeOld);
574 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
575 if (hlttreeOld) hltesd->ReadFromTree(hlttreeOld);
579 // create the ESD output file and tree
580 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
581 file->SetCompressionLevel(2);
582 if (!file->IsOpen()) {
583 AliError("opening AliESDs.root failed");
584 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
587 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
589 esd->CreateStdContent();
590 esd->WriteToTree(tree);
592 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
593 hltesd = new AliESD();
594 hltesd->CreateStdContent();
595 hltesd->WriteToTree(hlttree);
598 delete esd; delete hltesd;
599 esd = NULL; hltesd = NULL;
601 // create the branch with ESD additions
602 AliESDfriend *esdf = 0;
603 if (fWriteESDfriend) {
604 esdf = new AliESDfriend();
605 TBranch *br=tree->Branch("ESDfriend.","AliESDfriend", &esdf);
606 br->SetFile("AliESDfriends.root");
607 esd->AddObject(esdf);
610 // Get the diamond profile from OCDB
611 AliCDBEntry* entry = AliCDBManager::Instance()
612 ->Get("GRP/Calib/MeanVertex");
615 fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
617 AliError("No diamond profile found in OCDB!");
620 AliVertexerTracks tVertexer(AliTracker::GetBz());
621 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
624 if (fRawReader) fRawReader->RewindEvents();
626 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
627 if (fRawReader) fRawReader->NextEvent();
628 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
629 // copy old ESD to the new one
631 esd->ReadFromTree(treeOld);
632 treeOld->GetEntry(iEvent);
636 esd->ReadFromTree(hlttreeOld);
637 hlttreeOld->GetEntry(iEvent);
643 AliInfo(Form("processing event %d", iEvent));
644 fRunLoader->GetEvent(iEvent);
647 sprintf(aFileName, "ESD_%d.%d_final.root",
648 fRunLoader->GetHeader()->GetRun(),
649 fRunLoader->GetHeader()->GetEventNrInRun());
650 if (!gSystem->AccessPathName(aFileName)) continue;
652 // local reconstruction
653 if (!fRunLocalReconstruction.IsNull()) {
654 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
655 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
660 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
661 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
662 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
663 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
665 // Set magnetic field from the tracker
666 esd->SetMagneticField(AliTracker::GetBz());
667 hltesd->SetMagneticField(AliTracker::GetBz());
671 // Fill raw-data error log into the ESD
672 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
675 if (fRunVertexFinder) {
676 if (!ReadESD(esd, "vertex")) {
677 if (!RunVertexFinder(esd)) {
678 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
680 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
685 if (!fRunTracking.IsNull()) {
686 if (fRunHLTTracking) {
687 hltesd->SetVertex(esd->GetVertex());
688 if (!RunHLTTracking(hltesd)) {
689 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
695 if (!fRunTracking.IsNull()) {
696 if (fRunMuonTracking) {
697 if (!RunMuonTracking()) {
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;}
719 // fill Event header information from the RawEventHeader
720 if (fRawReader){FillRawEventHeaderESD(esd);}
723 AliESDpid::MakePID(esd);
724 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
726 if (fFillTriggerESD) {
727 if (!ReadESD(esd, "trigger")) {
728 if (!FillTriggerESD(esd)) {
729 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
731 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
735 //Try to improve the reconstructed primary vertex position using the tracks
736 AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
737 if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
740 if (pvtx->GetStatus()) {
741 // Store the improved primary vertex
742 esd->SetPrimaryVertex(pvtx);
743 // Propagate the tracks to the DCA to the improved primary vertex
744 Double_t somethingbig = 777.;
745 Double_t bz = esd->GetMagneticField();
746 Int_t nt=esd->GetNumberOfTracks();
748 AliESDtrack *t = esd->GetTrack(nt);
749 t->RelateToVertex(pvtx, bz, somethingbig);
756 vtxer.Tracks2V0vertices(esd);
759 AliCascadeVertexer cvtxer;
760 cvtxer.V0sTracks2CascadeVertices(esd);
764 if (fWriteESDfriend) {
765 new (esdf) AliESDfriend(); // Reset...
766 esd->GetESDfriend(esdf);
773 if (fCheckPointLevel > 0) WriteESD(esd, "final");
777 // delete esdf; esdf = 0;
781 tree->GetUserInfo()->Add(esd);
782 hlttree->GetUserInfo()->Add(hltesd);
786 if(fESDPar.Contains("ESD.par")){
787 AliInfo("Attaching ESD.par to Tree");
788 TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
789 tree->GetUserInfo()->Add(fn);
793 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
794 stopwatch.RealTime(),stopwatch.CpuTime()));
798 tree->SetBranchStatus("ESDfriend*",0);
803 TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
804 ESDFile2AODFile(file, aodFile);
809 // Create tags for the events in the ESD tree (the ESD tree is always present)
810 // In case of empty events the tags will contain dummy values
813 CleanUp(file, fileOld);
819 //_____________________________________________________________________________
820 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
822 // run the local reconstruction
824 TStopwatch stopwatch;
827 AliCDBManager* man = AliCDBManager::Instance();
828 Bool_t origCache = man->GetCacheFlag();
830 TString detStr = detectors;
831 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
832 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
833 AliReconstructor* reconstructor = GetReconstructor(iDet);
834 if (!reconstructor) continue;
835 if (reconstructor->HasLocalReconstruction()) continue;
837 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
838 TStopwatch stopwatchDet;
839 stopwatchDet.Start();
841 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
843 man->SetCacheFlag(kTRUE);
844 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
845 man->GetAll(calibPath); // entries are cached!
848 fRawReader->RewindEvents();
849 reconstructor->Reconstruct(fRunLoader, fRawReader);
851 reconstructor->Reconstruct(fRunLoader);
853 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
854 fgkDetectorName[iDet],
855 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
857 // unload calibration data
861 man->SetCacheFlag(origCache);
863 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
864 AliError(Form("the following detectors were not found: %s",
866 if (fStopOnError) return kFALSE;
869 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
870 stopwatch.RealTime(),stopwatch.CpuTime()));
875 //_____________________________________________________________________________
876 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
878 // run the local reconstruction
880 TStopwatch stopwatch;
883 TString detStr = detectors;
884 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
885 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
886 AliReconstructor* reconstructor = GetReconstructor(iDet);
887 if (!reconstructor) continue;
888 AliLoader* loader = fLoader[iDet];
890 // conversion of digits
891 if (fRawReader && reconstructor->HasDigitConversion()) {
892 AliInfo(Form("converting raw data digits into root objects for %s",
893 fgkDetectorName[iDet]));
894 TStopwatch stopwatchDet;
895 stopwatchDet.Start();
896 loader->LoadDigits("update");
897 loader->CleanDigits();
898 loader->MakeDigitsContainer();
899 TTree* digitsTree = loader->TreeD();
900 reconstructor->ConvertDigits(fRawReader, digitsTree);
901 loader->WriteDigits("OVERWRITE");
902 loader->UnloadDigits();
903 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
904 fgkDetectorName[iDet],
905 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
908 // local reconstruction
909 if (!reconstructor->HasLocalReconstruction()) continue;
910 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
911 TStopwatch stopwatchDet;
912 stopwatchDet.Start();
913 loader->LoadRecPoints("update");
914 loader->CleanRecPoints();
915 loader->MakeRecPointsContainer();
916 TTree* clustersTree = loader->TreeR();
917 if (fRawReader && !reconstructor->HasDigitConversion()) {
918 reconstructor->Reconstruct(fRawReader, clustersTree);
920 loader->LoadDigits("read");
921 TTree* digitsTree = loader->TreeD();
923 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
924 if (fStopOnError) return kFALSE;
926 reconstructor->Reconstruct(digitsTree, clustersTree);
928 loader->UnloadDigits();
930 loader->WriteRecPoints("OVERWRITE");
931 loader->UnloadRecPoints();
932 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
933 fgkDetectorName[iDet],
934 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
937 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
938 AliError(Form("the following detectors were not found: %s",
940 if (fStopOnError) return kFALSE;
943 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
944 stopwatch.RealTime(),stopwatch.CpuTime()));
949 //_____________________________________________________________________________
950 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
952 // run the barrel tracking
954 TStopwatch stopwatch;
957 AliESDVertex* vertex = NULL;
958 Double_t vtxPos[3] = {0, 0, 0};
959 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
961 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
962 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
963 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
967 if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
968 AliInfo("running the ITS vertex finder");
969 if (fLoader[0]) fLoader[0]->LoadRecPoints();
970 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
971 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
973 AliWarning("Vertex not found");
974 vertex = new AliESDVertex();
975 vertex->SetName("default");
978 vertex->SetName("reconstructed");
982 AliInfo("getting the primary vertex from MC");
983 vertex = new AliESDVertex(vtxPos, vtxErr);
987 vertex->GetXYZ(vtxPos);
988 vertex->GetSigmaXYZ(vtxErr);
990 AliWarning("no vertex reconstructed");
991 vertex = new AliESDVertex(vtxPos, vtxErr);
993 esd->SetVertex(vertex);
994 // if SPD multiplicity has been determined, it is stored in the ESD
995 AliMultiplicity *mult = fVertexer->GetMultiplicity();
996 if(mult)esd->SetMultiplicity(mult);
998 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
999 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1003 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1004 stopwatch.RealTime(),stopwatch.CpuTime()));
1009 //_____________________________________________________________________________
1010 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1012 // run the HLT barrel tracking
1014 TStopwatch stopwatch;
1018 AliError("Missing runLoader!");
1022 AliInfo("running HLT tracking");
1024 // Get a pointer to the HLT reconstructor
1025 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1026 if (!reconstructor) return kFALSE;
1029 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1030 TString detName = fgkDetectorName[iDet];
1031 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1032 reconstructor->SetOption(detName.Data());
1033 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1035 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1036 if (fStopOnError) return kFALSE;
1040 Double_t vtxErr[3]={0.005,0.005,0.010};
1041 const AliESDVertex *vertex = esd->GetVertex();
1042 vertex->GetXYZ(vtxPos);
1043 tracker->SetVertex(vtxPos,vtxErr);
1045 fLoader[iDet]->LoadRecPoints("read");
1046 TTree* tree = fLoader[iDet]->TreeR();
1048 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1051 tracker->LoadClusters(tree);
1053 if (tracker->Clusters2Tracks(esd) != 0) {
1054 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1058 tracker->UnloadClusters();
1063 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1064 stopwatch.RealTime(),stopwatch.CpuTime()));
1069 //_____________________________________________________________________________
1070 Bool_t AliReconstruction::RunMuonTracking()
1072 // run the muon spectrometer tracking
1074 TStopwatch stopwatch;
1078 AliError("Missing runLoader!");
1081 Int_t iDet = 7; // for MUON
1083 AliInfo("is running...");
1085 // Get a pointer to the MUON reconstructor
1086 AliReconstructor *reconstructor = GetReconstructor(iDet);
1087 if (!reconstructor) return kFALSE;
1090 TString detName = fgkDetectorName[iDet];
1091 AliDebug(1, Form("%s tracking", detName.Data()));
1092 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1094 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1099 fLoader[iDet]->LoadTracks("update");
1100 fLoader[iDet]->CleanTracks();
1101 fLoader[iDet]->MakeTracksContainer();
1104 fLoader[iDet]->LoadRecPoints("read");
1106 if (!tracker->Clusters2Tracks(0x0)) {
1107 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1110 fLoader[iDet]->UnloadRecPoints();
1112 fLoader[iDet]->WriteTracks("OVERWRITE");
1113 fLoader[iDet]->UnloadTracks();
1118 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1119 stopwatch.RealTime(),stopwatch.CpuTime()));
1125 //_____________________________________________________________________________
1126 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1128 // run the barrel tracking
1130 TStopwatch stopwatch;
1133 AliInfo("running tracking");
1135 //Fill the ESD with the T0 info (will be used by the TOF)
1136 if (fReconstructor[11])
1137 GetReconstructor(11)->FillESD(fRunLoader, esd);
1139 // pass 1: TPC + ITS inwards
1140 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1141 if (!fTracker[iDet]) continue;
1142 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1145 fLoader[iDet]->LoadRecPoints("read");
1146 TTree* tree = fLoader[iDet]->TreeR();
1148 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1151 fTracker[iDet]->LoadClusters(tree);
1154 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1155 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1158 if (fCheckPointLevel > 1) {
1159 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1161 // preliminary PID in TPC needed by the ITS tracker
1163 GetReconstructor(1)->FillESD(fRunLoader, esd);
1164 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1165 AliESDpid::MakePID(esd);
1169 // pass 2: ALL backwards
1170 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1171 if (!fTracker[iDet]) continue;
1172 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1175 if (iDet > 1) { // all except ITS, TPC
1177 fLoader[iDet]->LoadRecPoints("read");
1178 tree = fLoader[iDet]->TreeR();
1180 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1183 fTracker[iDet]->LoadClusters(tree);
1187 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1188 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1191 if (fCheckPointLevel > 1) {
1192 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1196 if (iDet > 2) { // all except ITS, TPC, TRD
1197 fTracker[iDet]->UnloadClusters();
1198 fLoader[iDet]->UnloadRecPoints();
1200 // updated PID in TPC needed by the ITS tracker -MI
1202 GetReconstructor(1)->FillESD(fRunLoader, esd);
1203 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1204 AliESDpid::MakePID(esd);
1208 // write space-points to the ESD in case alignment data output
1210 if (fWriteAlignmentData)
1211 WriteAlignmentData(esd);
1213 // pass 3: TRD + TPC + ITS refit inwards
1214 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1215 if (!fTracker[iDet]) continue;
1216 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1219 if (fTracker[iDet]->RefitInward(esd) != 0) {
1220 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1223 if (fCheckPointLevel > 1) {
1224 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1228 fTracker[iDet]->UnloadClusters();
1229 fLoader[iDet]->UnloadRecPoints();
1232 // Propagate track to the vertex - if not done by ITS
1234 Int_t ntracks = esd->GetNumberOfTracks();
1235 for (Int_t itrack=0; itrack<ntracks; itrack++){
1236 const Double_t kRadius = 3; // beam pipe radius
1237 const Double_t kMaxStep = 5; // max step
1238 const Double_t kMaxD = 123456; // max distance to prim vertex
1239 Double_t fieldZ = AliTracker::GetBz(); //
1240 AliESDtrack * track = esd->GetTrack(itrack);
1241 if (!track) continue;
1242 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1243 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1244 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1247 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1248 stopwatch.RealTime(),stopwatch.CpuTime()));
1253 //_____________________________________________________________________________
1254 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1256 // fill the event summary data
1258 TStopwatch stopwatch;
1260 AliInfo("filling ESD");
1262 TString detStr = detectors;
1263 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1264 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1265 AliReconstructor* reconstructor = GetReconstructor(iDet);
1266 if (!reconstructor) continue;
1268 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1269 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1270 TTree* clustersTree = NULL;
1271 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1272 fLoader[iDet]->LoadRecPoints("read");
1273 clustersTree = fLoader[iDet]->TreeR();
1274 if (!clustersTree) {
1275 AliError(Form("Can't get the %s clusters tree",
1276 fgkDetectorName[iDet]));
1277 if (fStopOnError) return kFALSE;
1280 if (fRawReader && !reconstructor->HasDigitConversion()) {
1281 reconstructor->FillESD(fRawReader, clustersTree, esd);
1283 TTree* digitsTree = NULL;
1284 if (fLoader[iDet]) {
1285 fLoader[iDet]->LoadDigits("read");
1286 digitsTree = fLoader[iDet]->TreeD();
1288 AliError(Form("Can't get the %s digits tree",
1289 fgkDetectorName[iDet]));
1290 if (fStopOnError) return kFALSE;
1293 reconstructor->FillESD(digitsTree, clustersTree, esd);
1294 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1296 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1297 fLoader[iDet]->UnloadRecPoints();
1301 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1303 reconstructor->FillESD(fRunLoader, esd);
1305 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1309 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1310 AliError(Form("the following detectors were not found: %s",
1312 if (fStopOnError) return kFALSE;
1315 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1316 stopwatch.RealTime(),stopwatch.CpuTime()));
1321 //_____________________________________________________________________________
1322 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1324 // Reads the trigger decision which is
1325 // stored in Trigger.root file and fills
1326 // the corresponding esd entries
1328 AliInfo("Filling trigger information into the ESD");
1331 AliCTPRawStream input(fRawReader);
1332 if (!input.Next()) {
1333 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1336 esd->SetTriggerMask(input.GetClassMask());
1337 esd->SetTriggerCluster(input.GetClusterMask());
1340 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1342 if (!runloader->LoadTrigger()) {
1343 AliCentralTrigger *aCTP = runloader->GetTrigger();
1344 esd->SetTriggerMask(aCTP->GetClassMask());
1345 esd->SetTriggerCluster(aCTP->GetClusterMask());
1348 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1353 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1365 //_____________________________________________________________________________
1366 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1369 // Filling information from RawReader Header
1372 AliInfo("Filling information from RawReader Header");
1373 esd->SetBunchCrossNumber(0);
1374 esd->SetOrbitNumber(0);
1375 esd->SetPeriodNumber(0);
1376 esd->SetTimeStamp(0);
1377 esd->SetEventType(0);
1378 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1381 const UInt_t *id = eventHeader->GetP("Id");
1382 esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1383 esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1384 esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1386 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1387 esd->SetEventType((eventHeader->Get("Type")));
1394 //_____________________________________________________________________________
1395 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1397 // check whether detName is contained in detectors
1398 // if yes, it is removed from detectors
1400 // check if all detectors are selected
1401 if ((detectors.CompareTo("ALL") == 0) ||
1402 detectors.BeginsWith("ALL ") ||
1403 detectors.EndsWith(" ALL") ||
1404 detectors.Contains(" ALL ")) {
1409 // search for the given detector
1410 Bool_t result = kFALSE;
1411 if ((detectors.CompareTo(detName) == 0) ||
1412 detectors.BeginsWith(detName+" ") ||
1413 detectors.EndsWith(" "+detName) ||
1414 detectors.Contains(" "+detName+" ")) {
1415 detectors.ReplaceAll(detName, "");
1419 // clean up the detectors string
1420 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1421 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1422 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1427 //_____________________________________________________________________________
1428 Bool_t AliReconstruction::InitRunLoader()
1430 // get or create the run loader
1432 if (gAlice) delete gAlice;
1435 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1436 // load all base libraries to get the loader classes
1437 TString libs = gSystem->GetLibraries();
1438 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1439 TString detName = fgkDetectorName[iDet];
1440 if (detName == "HLT") continue;
1441 if (libs.Contains("lib" + detName + "base.so")) continue;
1442 gSystem->Load("lib" + detName + "base.so");
1444 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1446 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1450 fRunLoader->CdGAFile();
1451 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1452 if (fRunLoader->LoadgAlice() == 0) {
1453 gAlice = fRunLoader->GetAliRun();
1454 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1457 if (!gAlice && !fRawReader) {
1458 AliError(Form("no gAlice object found in file %s",
1459 fGAliceFileName.Data()));
1464 } else { // galice.root does not exist
1466 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1470 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1471 AliConfig::GetDefaultEventFolderName(),
1474 AliError(Form("could not create run loader in file %s",
1475 fGAliceFileName.Data()));
1479 fRunLoader->MakeTree("E");
1481 while (fRawReader->NextEvent()) {
1482 fRunLoader->SetEventNumber(iEvent);
1483 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1485 fRunLoader->MakeTree("H");
1486 fRunLoader->TreeE()->Fill();
1489 fRawReader->RewindEvents();
1490 if (fNumberOfEventsPerFile > 0)
1491 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1493 fRunLoader->SetNumberOfEventsPerFile(iEvent);
1494 fRunLoader->WriteHeader("OVERWRITE");
1495 fRunLoader->CdGAFile();
1496 fRunLoader->Write(0, TObject::kOverwrite);
1497 // AliTracker::SetFieldMap(???);
1503 //_____________________________________________________________________________
1504 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1506 // get the reconstructor object and the loader for a detector
1508 if (fReconstructor[iDet]) return fReconstructor[iDet];
1510 // load the reconstructor object
1511 TPluginManager* pluginManager = gROOT->GetPluginManager();
1512 TString detName = fgkDetectorName[iDet];
1513 TString recName = "Ali" + detName + "Reconstructor";
1514 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1516 if (detName == "HLT") {
1517 if (!gROOT->GetClass("AliLevel3")) {
1518 gSystem->Load("libAliHLTSrc.so");
1519 gSystem->Load("libAliHLTMisc.so");
1520 gSystem->Load("libAliHLTHough.so");
1521 gSystem->Load("libAliHLTComp.so");
1525 AliReconstructor* reconstructor = NULL;
1526 // first check if a plugin is defined for the reconstructor
1527 TPluginHandler* pluginHandler =
1528 pluginManager->FindHandler("AliReconstructor", detName);
1529 // if not, add a plugin for it
1530 if (!pluginHandler) {
1531 AliDebug(1, Form("defining plugin for %s", recName.Data()));
1532 TString libs = gSystem->GetLibraries();
1533 if (libs.Contains("lib" + detName + "base.so") ||
1534 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1535 pluginManager->AddHandler("AliReconstructor", detName,
1536 recName, detName + "rec", recName + "()");
1538 pluginManager->AddHandler("AliReconstructor", detName,
1539 recName, detName, recName + "()");
1541 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1543 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1544 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1546 if (reconstructor) {
1547 TObject* obj = fOptions.FindObject(detName.Data());
1548 if (obj) reconstructor->SetOption(obj->GetTitle());
1549 reconstructor->Init(fRunLoader);
1550 fReconstructor[iDet] = reconstructor;
1553 // get or create the loader
1554 if (detName != "HLT") {
1555 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1556 if (!fLoader[iDet]) {
1557 AliConfig::Instance()
1558 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1560 // first check if a plugin is defined for the loader
1562 pluginManager->FindHandler("AliLoader", detName);
1563 // if not, add a plugin for it
1564 if (!pluginHandler) {
1565 TString loaderName = "Ali" + detName + "Loader";
1566 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1567 pluginManager->AddHandler("AliLoader", detName,
1568 loaderName, detName + "base",
1569 loaderName + "(const char*, TFolder*)");
1570 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1572 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1574 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1575 fRunLoader->GetEventFolder());
1577 if (!fLoader[iDet]) { // use default loader
1578 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1580 if (!fLoader[iDet]) {
1581 AliWarning(Form("couldn't get loader for %s", detName.Data()));
1582 if (fStopOnError) return NULL;
1584 fRunLoader->AddLoader(fLoader[iDet]);
1585 fRunLoader->CdGAFile();
1586 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1587 fRunLoader->Write(0, TObject::kOverwrite);
1592 return reconstructor;
1595 //_____________________________________________________________________________
1596 Bool_t AliReconstruction::CreateVertexer()
1598 // create the vertexer
1601 AliReconstructor* itsReconstructor = GetReconstructor(0);
1602 if (itsReconstructor) {
1603 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1606 AliWarning("couldn't create a vertexer for ITS");
1607 if (fStopOnError) return kFALSE;
1613 //_____________________________________________________________________________
1614 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1616 // create the trackers
1618 TString detStr = detectors;
1619 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1620 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1621 AliReconstructor* reconstructor = GetReconstructor(iDet);
1622 if (!reconstructor) continue;
1623 TString detName = fgkDetectorName[iDet];
1624 if (detName == "HLT") {
1625 fRunHLTTracking = kTRUE;
1628 if (detName == "MUON") {
1629 fRunMuonTracking = kTRUE;
1634 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1635 if (!fTracker[iDet] && (iDet < 7)) {
1636 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1637 if (fStopOnError) return kFALSE;
1644 //_____________________________________________________________________________
1645 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1647 // delete trackers and the run loader and close and delete the file
1649 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1650 delete fReconstructor[iDet];
1651 fReconstructor[iDet] = NULL;
1652 fLoader[iDet] = NULL;
1653 delete fTracker[iDet];
1654 fTracker[iDet] = NULL;
1658 delete fDiamondProfile;
1659 fDiamondProfile = NULL;
1674 gSystem->Unlink("AliESDs.old.root");
1679 //_____________________________________________________________________________
1680 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1682 // read the ESD event from a file
1684 if (!esd) return kFALSE;
1686 sprintf(fileName, "ESD_%d.%d_%s.root",
1687 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1688 if (gSystem->AccessPathName(fileName)) return kFALSE;
1690 AliInfo(Form("reading ESD from file %s", fileName));
1691 AliDebug(1, Form("reading ESD from file %s", fileName));
1692 TFile* file = TFile::Open(fileName);
1693 if (!file || !file->IsOpen()) {
1694 AliError(Form("opening %s failed", fileName));
1701 esd = (AliESD*) file->Get("ESD");
1707 //_____________________________________________________________________________
1708 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1710 // write the ESD event to a file
1714 sprintf(fileName, "ESD_%d.%d_%s.root",
1715 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1717 AliDebug(1, Form("writing ESD to file %s", fileName));
1718 TFile* file = TFile::Open(fileName, "recreate");
1719 if (!file || !file->IsOpen()) {
1720 AliError(Form("opening %s failed", fileName));
1731 //_____________________________________________________________________________
1732 void AliReconstruction::CreateTag(TFile* file)
1735 Float_t lhcLuminosity = 0.0;
1736 TString lhcState = "test";
1737 UInt_t detectorMask = 0;
1742 Double_t fMUONMASS = 0.105658369;
1745 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1746 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1748 TLorentzVector fEPvector;
1750 Float_t fZVertexCut = 10.0;
1751 Float_t fRhoVertexCut = 2.0;
1753 Float_t fLowPtCut = 1.0;
1754 Float_t fHighPtCut = 3.0;
1755 Float_t fVeryHighPtCut = 10.0;
1758 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1760 // Creates the tags for all the events in a given ESD file
1762 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1763 Int_t nPos, nNeg, nNeutr;
1764 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1765 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1766 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1767 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1768 Float_t maxPt = .0, meanPt = .0, totalP = .0;
1770 Int_t iRunNumber = 0;
1771 TString fVertexName("default");
1773 AliRunTag *tag = new AliRunTag();
1774 AliEventTag *evTag = new AliEventTag();
1775 TTree ttag("T","A Tree with event tags");
1776 TBranch * btag = ttag.Branch("AliTAG", &tag);
1777 btag->SetCompressionLevel(9);
1779 AliInfo(Form("Creating the tags......."));
1781 if (!file || !file->IsOpen()) {
1782 AliError(Form("opening failed"));
1786 Int_t lastEvent = 0;
1787 TTree *b = (TTree*) file->Get("esdTree");
1788 AliESD *esd = new AliESD();
1789 esd->ReadFromTree(b);
1791 b->GetEntry(fFirstEvent);
1792 Int_t iInitRunNumber = esd->GetRunNumber();
1794 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1795 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1796 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1824 b->GetEntry(iEventNumber);
1825 iRunNumber = esd->GetRunNumber();
1826 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1827 const AliESDVertex * vertexIn = esd->GetVertex();
1828 if (!vertexIn) AliError("ESD has not defined vertex.");
1829 if (vertexIn) fVertexName = vertexIn->GetName();
1830 if(fVertexName != "default") fVertexflag = 1;
1831 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1832 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1833 UInt_t status = esdTrack->GetStatus();
1835 //select only tracks with ITS refit
1836 if ((status&AliESDtrack::kITSrefit)==0) continue;
1837 //select only tracks with TPC refit
1838 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1840 //select only tracks with the "combined PID"
1841 if ((status&AliESDtrack::kESDpid)==0) continue;
1843 esdTrack->GetPxPyPz(p);
1844 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1845 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1848 if(fPt > maxPt) maxPt = fPt;
1850 if(esdTrack->GetSign() > 0) {
1852 if(fPt > fLowPtCut) nCh1GeV++;
1853 if(fPt > fHighPtCut) nCh3GeV++;
1854 if(fPt > fVeryHighPtCut) nCh10GeV++;
1856 if(esdTrack->GetSign() < 0) {
1858 if(fPt > fLowPtCut) nCh1GeV++;
1859 if(fPt > fHighPtCut) nCh3GeV++;
1860 if(fPt > fVeryHighPtCut) nCh10GeV++;
1862 if(esdTrack->GetSign() == 0) nNeutr++;
1866 esdTrack->GetESDpid(prob);
1869 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1870 if(rcc == 0.0) continue;
1873 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1876 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1878 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1880 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1882 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1884 if(fPt > fLowPtCut) nEl1GeV++;
1885 if(fPt > fHighPtCut) nEl3GeV++;
1886 if(fPt > fVeryHighPtCut) nEl10GeV++;
1894 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1895 // loop over all reconstructed tracks (also first track of combination)
1896 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1897 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1898 if (muonTrack == 0x0) continue;
1900 // Coordinates at vertex
1901 fZ = muonTrack->GetZ();
1902 fY = muonTrack->GetBendingCoor();
1903 fX = muonTrack->GetNonBendingCoor();
1905 fThetaX = muonTrack->GetThetaX();
1906 fThetaY = muonTrack->GetThetaY();
1908 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1909 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1910 fPxRec = fPzRec * TMath::Tan(fThetaX);
1911 fPyRec = fPzRec * TMath::Tan(fThetaY);
1912 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1914 //ChiSquare of the track if needed
1915 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1916 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1917 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1919 // total number of muons inside a vertex cut
1920 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1922 if(fEPvector.Pt() > fLowPtCut) {
1924 if(fEPvector.Pt() > fHighPtCut) {
1926 if (fEPvector.Pt() > fVeryHighPtCut) {
1934 // Fill the event tags
1936 meanPt = meanPt/ntrack;
1938 evTag->SetEventId(iEventNumber+1);
1940 evTag->SetVertexX(vertexIn->GetXv());
1941 evTag->SetVertexY(vertexIn->GetYv());
1942 evTag->SetVertexZ(vertexIn->GetZv());
1943 evTag->SetVertexZError(vertexIn->GetZRes());
1945 evTag->SetVertexFlag(fVertexflag);
1947 evTag->SetT0VertexZ(esd->GetT0zVertex());
1949 evTag->SetTriggerMask(esd->GetTriggerMask());
1950 evTag->SetTriggerCluster(esd->GetTriggerCluster());
1952 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1953 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1954 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1955 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1956 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1957 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1960 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1961 evTag->SetNumOfPosTracks(nPos);
1962 evTag->SetNumOfNegTracks(nNeg);
1963 evTag->SetNumOfNeutrTracks(nNeutr);
1965 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1966 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1967 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1968 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1970 evTag->SetNumOfProtons(nProtons);
1971 evTag->SetNumOfKaons(nKaons);
1972 evTag->SetNumOfPions(nPions);
1973 evTag->SetNumOfMuons(nMuons);
1974 evTag->SetNumOfElectrons(nElectrons);
1975 evTag->SetNumOfPhotons(nGamas);
1976 evTag->SetNumOfPi0s(nPi0s);
1977 evTag->SetNumOfNeutrons(nNeutrons);
1978 evTag->SetNumOfKaon0s(nK0s);
1980 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1981 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1982 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1983 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1984 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1985 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1986 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1987 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1988 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1990 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1991 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1993 evTag->SetTotalMomentum(totalP);
1994 evTag->SetMeanPt(meanPt);
1995 evTag->SetMaxPt(maxPt);
1997 tag->SetLHCTag(lhcLuminosity,lhcState);
1998 tag->SetDetectorTag(detectorMask);
2000 tag->SetRunId(iInitRunNumber);
2001 tag->AddEventTag(*evTag);
2003 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2004 else lastEvent = fLastEvent;
2010 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2011 tag->GetRunId(),fFirstEvent,lastEvent );
2012 AliInfo(Form("writing tags to file %s", fileName));
2013 AliDebug(1, Form("writing tags to file %s", fileName));
2015 TFile* ftag = TFile::Open(fileName, "recreate");
2024 //_____________________________________________________________________________
2025 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
2027 // write all files from the given esd file to an aod file
2029 // create an AliAOD object
2030 AliAODEvent *aod = new AliAODEvent();
2031 aod->CreateStdContent();
2037 TTree *aodTree = new TTree("AOD", "AliAOD tree");
2038 aodTree->Branch(aod->GetList());
2041 TTree *t = (TTree*) esdFile->Get("esdTree");
2042 TBranch *b = t->GetBranch("ESD");
2044 b->SetAddress(&esd);
2046 Int_t nEvents = b->GetEntries();
2048 // set arrays and pointers
2056 // loop over events and fill them
2057 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2058 b->GetEntry(iEvent);
2060 // Multiplicity information needed by the header (to be revised!)
2061 Int_t nTracks = esd->GetNumberOfTracks();
2062 Int_t nPosTracks = 0;
2063 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2064 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2066 // create the header
2067 aod->AddHeader(new AliAODHeader(esd->GetRunNumber(),
2068 esd->GetBunchCrossNumber(),
2069 esd->GetOrbitNumber(),
2070 esd->GetPeriodNumber(),
2074 esd->GetMagneticField(),
2075 -999., // fill muon magnetic field
2076 -999., // centrality; to be filled, still
2077 esd->GetZDCN1Energy(),
2078 esd->GetZDCP1Energy(),
2079 esd->GetZDCN2Energy(),
2080 esd->GetZDCP2Energy(),
2081 esd->GetZDCEMEnergy(),
2082 esd->GetTriggerMask(),
2083 esd->GetTriggerCluster(),
2084 esd->GetEventType()));
2086 Int_t nV0s = esd->GetNumberOfV0s();
2087 Int_t nCascades = esd->GetNumberOfCascades();
2088 Int_t nKinks = esd->GetNumberOfKinks();
2089 Int_t nVertices = nV0s + nCascades + nKinks;
2091 aod->ResetStd(nTracks, nVertices);
2092 AliAODTrack *aodTrack;
2095 // Array to take into account the tracks already added to the AOD
2096 Bool_t * usedTrack = NULL;
2098 usedTrack = new Bool_t[nTracks];
2099 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2101 // Array to take into account the V0s already added to the AOD
2102 Bool_t * usedV0 = NULL;
2104 usedV0 = new Bool_t[nV0s];
2105 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2107 // Array to take into account the kinks already added to the AOD
2108 Bool_t * usedKink = NULL;
2110 usedKink = new Bool_t[nKinks];
2111 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2114 // Access to the AOD container of vertices
2115 TClonesArray &vertices = *(aod->GetVertices());
2118 // Access to the AOD container of tracks
2119 TClonesArray &tracks = *(aod->GetTracks());
2122 // Add primary vertex. The primary tracks will be defined
2123 // after the loops on the composite objects (V0, cascades, kinks)
2124 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2126 vtx->GetXYZ(pos); // position
2127 vtx->GetCovMatrix(covVtx); //covariance matrix
2129 AliAODVertex * primary = new(vertices[jVertices++])
2130 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2132 // Create vertices starting from the most complex objects
2135 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2136 AliESDcascade *cascade = esd->GetCascade(nCascade);
2138 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2139 cascade->GetPosCovXi(covVtx);
2141 // Add the cascade vertex
2142 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2144 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2146 AliAODVertex::kCascade);
2148 primary->AddDaughter(vcascade);
2150 // Add the V0 from the cascade. The ESD class have to be optimized...
2151 // Now we have to search for the corresponding Vo in the list of V0s
2152 // using the indeces of the positive and negative tracks
2154 Int_t posFromV0 = cascade->GetPindex();
2155 Int_t negFromV0 = cascade->GetNindex();
2158 AliESDv0 * v0 = 0x0;
2161 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2163 v0 = esd->GetV0(iV0);
2164 Int_t posV0 = v0->GetPindex();
2165 Int_t negV0 = v0->GetNindex();
2167 if (posV0==posFromV0 && negV0==negFromV0) {
2173 AliAODVertex * vV0FromCascade = 0x0;
2175 if (indV0>-1 && !usedV0[indV0] ) {
2177 // the V0 exists in the array of V0s and is not used
2179 usedV0[indV0] = kTRUE;
2181 v0->GetXYZ(pos[0], pos[1], pos[2]);
2182 v0->GetPosCov(covVtx);
2184 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2186 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2191 // the V0 doesn't exist in the array of V0s or was used
2192 cerr << "Error: event " << iEvent << " cascade " << nCascade
2193 << " The V0 " << indV0
2194 << " doesn't exist in the array of V0s or was used!" << endl;
2196 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2197 cascade->GetPosCov(covVtx);
2199 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2201 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2204 vcascade->AddDaughter(vV0FromCascade);
2207 // Add the positive tracks from the V0
2209 if (! usedTrack[posFromV0]) {
2211 usedTrack[posFromV0] = kTRUE;
2213 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2214 esdTrack->GetPxPyPz(p);
2215 esdTrack->GetXYZ(pos);
2216 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2217 esdTrack->GetESDpid(pid);
2219 vV0FromCascade->AddDaughter(aodTrack =
2220 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2221 esdTrack->GetLabel(),
2227 (Short_t)esdTrack->GetSign(),
2228 esdTrack->GetITSClusterMap(),
2231 kTRUE, // check if this is right
2232 kFALSE, // check if this is right
2233 AliAODTrack::kSecondary)
2235 aodTrack->ConvertAliPIDtoAODPID();
2238 cerr << "Error: event " << iEvent << " cascade " << nCascade
2239 << " track " << posFromV0 << " has already been used!" << endl;
2242 // Add the negative tracks from the V0
2244 if (!usedTrack[negFromV0]) {
2246 usedTrack[negFromV0] = kTRUE;
2248 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2249 esdTrack->GetPxPyPz(p);
2250 esdTrack->GetXYZ(pos);
2251 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2252 esdTrack->GetESDpid(pid);
2254 vV0FromCascade->AddDaughter(aodTrack =
2255 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2256 esdTrack->GetLabel(),
2262 (Short_t)esdTrack->GetSign(),
2263 esdTrack->GetITSClusterMap(),
2266 kTRUE, // check if this is right
2267 kFALSE, // check if this is right
2268 AliAODTrack::kSecondary)
2270 aodTrack->ConvertAliPIDtoAODPID();
2273 cerr << "Error: event " << iEvent << " cascade " << nCascade
2274 << " track " << negFromV0 << " has already been used!" << endl;
2277 // Add the bachelor track from the cascade
2279 Int_t bachelor = cascade->GetBindex();
2281 if(!usedTrack[bachelor]) {
2283 usedTrack[bachelor] = kTRUE;
2285 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2286 esdTrack->GetPxPyPz(p);
2287 esdTrack->GetXYZ(pos);
2288 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2289 esdTrack->GetESDpid(pid);
2291 vcascade->AddDaughter(aodTrack =
2292 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2293 esdTrack->GetLabel(),
2299 (Short_t)esdTrack->GetSign(),
2300 esdTrack->GetITSClusterMap(),
2303 kTRUE, // check if this is right
2304 kFALSE, // check if this is right
2305 AliAODTrack::kSecondary)
2307 aodTrack->ConvertAliPIDtoAODPID();
2310 cerr << "Error: event " << iEvent << " cascade " << nCascade
2311 << " track " << bachelor << " has already been used!" << endl;
2314 // Add the primary track of the cascade (if any)
2316 } // end of the loop on cascades
2320 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2322 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2324 AliESDv0 *v0 = esd->GetV0(nV0);
2326 v0->GetXYZ(pos[0], pos[1], pos[2]);
2327 v0->GetPosCov(covVtx);
2329 AliAODVertex * vV0 =
2330 new(vertices[jVertices++]) AliAODVertex(pos,
2332 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2335 primary->AddDaughter(vV0);
2337 Int_t posFromV0 = v0->GetPindex();
2338 Int_t negFromV0 = v0->GetNindex();
2340 // Add the positive tracks from the V0
2342 if (!usedTrack[posFromV0]) {
2344 usedTrack[posFromV0] = kTRUE;
2346 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2347 esdTrack->GetPxPyPz(p);
2348 esdTrack->GetXYZ(pos);
2349 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2350 esdTrack->GetESDpid(pid);
2352 vV0->AddDaughter(aodTrack =
2353 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2354 esdTrack->GetLabel(),
2360 (Short_t)esdTrack->GetSign(),
2361 esdTrack->GetITSClusterMap(),
2364 kTRUE, // check if this is right
2365 kFALSE, // check if this is right
2366 AliAODTrack::kSecondary)
2368 aodTrack->ConvertAliPIDtoAODPID();
2371 cerr << "Error: event " << iEvent << " V0 " << nV0
2372 << " track " << posFromV0 << " has already been used!" << endl;
2375 // Add the negative tracks from the V0
2377 if (!usedTrack[negFromV0]) {
2379 usedTrack[negFromV0] = kTRUE;
2381 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2382 esdTrack->GetPxPyPz(p);
2383 esdTrack->GetXYZ(pos);
2384 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2385 esdTrack->GetESDpid(pid);
2387 vV0->AddDaughter(aodTrack =
2388 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2389 esdTrack->GetLabel(),
2395 (Short_t)esdTrack->GetSign(),
2396 esdTrack->GetITSClusterMap(),
2399 kTRUE, // check if this is right
2400 kFALSE, // check if this is right
2401 AliAODTrack::kSecondary)
2403 aodTrack->ConvertAliPIDtoAODPID();
2406 cerr << "Error: event " << iEvent << " V0 " << nV0
2407 << " track " << negFromV0 << " has already been used!" << endl;
2410 } // end of the loop on V0s
2412 // Kinks: it is a big mess the access to the information in the kinks
2413 // The loop is on the tracks in order to find the mother and daugther of each kink
2416 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2419 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2421 Int_t ikink = esdTrack->GetKinkIndex(0);
2424 // Negative kink index: mother, positive: daughter
2426 // Search for the second track of the kink
2428 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2430 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2432 Int_t jkink = esdTrack1->GetKinkIndex(0);
2434 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2436 // The two tracks are from the same kink
2438 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2441 Int_t idaughter = -1;
2443 if (ikink<0 && jkink>0) {
2448 else if (ikink>0 && jkink<0) {
2454 cerr << "Error: Wrong combination of kink indexes: "
2455 << ikink << " " << jkink << endl;
2459 // Add the mother track
2461 AliAODTrack * mother = NULL;
2463 if (!usedTrack[imother]) {
2465 usedTrack[imother] = kTRUE;
2467 AliESDtrack *esdTrack = esd->GetTrack(imother);
2468 esdTrack->GetPxPyPz(p);
2469 esdTrack->GetXYZ(pos);
2470 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2471 esdTrack->GetESDpid(pid);
2474 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2475 esdTrack->GetLabel(),
2481 (Short_t)esdTrack->GetSign(),
2482 esdTrack->GetITSClusterMap(),
2485 kTRUE, // check if this is right
2486 kTRUE, // check if this is right
2487 AliAODTrack::kPrimary);
2488 primary->AddDaughter(mother);
2489 mother->ConvertAliPIDtoAODPID();
2492 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2493 << " track " << imother << " has already been used!" << endl;
2496 // Add the kink vertex
2497 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2499 AliAODVertex * vkink =
2500 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2504 AliAODVertex::kKink);
2505 // Add the daughter track
2507 AliAODTrack * daughter = NULL;
2509 if (!usedTrack[idaughter]) {
2511 usedTrack[idaughter] = kTRUE;
2513 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2514 esdTrack->GetPxPyPz(p);
2515 esdTrack->GetXYZ(pos);
2516 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2517 esdTrack->GetESDpid(pid);
2520 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2521 esdTrack->GetLabel(),
2527 (Short_t)esdTrack->GetSign(),
2528 esdTrack->GetITSClusterMap(),
2531 kTRUE, // check if this is right
2532 kTRUE, // check if this is right
2533 AliAODTrack::kPrimary);
2534 vkink->AddDaughter(daughter);
2535 daughter->ConvertAliPIDtoAODPID();
2538 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2539 << " track " << idaughter << " has already been used!" << endl;
2551 // Tracks (primary and orphan)
2553 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2556 if (usedTrack[nTrack]) continue;
2558 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2559 esdTrack->GetPxPyPz(p);
2560 esdTrack->GetXYZ(pos);
2561 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2562 esdTrack->GetESDpid(pid);
2564 Float_t impactXY, impactZ;
2566 esdTrack->GetImpactParameters(impactXY,impactZ);
2569 // track inside the beam pipe
2571 primary->AddDaughter(aodTrack =
2572 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2573 esdTrack->GetLabel(),
2579 (Short_t)esdTrack->GetSign(),
2580 esdTrack->GetITSClusterMap(),
2583 kTRUE, // check if this is right
2584 kTRUE, // check if this is right
2585 AliAODTrack::kPrimary)
2587 aodTrack->ConvertAliPIDtoAODPID();
2590 // outside the beam pipe: orphan track
2592 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2593 esdTrack->GetLabel(),
2599 (Short_t)esdTrack->GetSign(),
2600 esdTrack->GetITSClusterMap(),
2603 kFALSE, // check if this is right
2604 kFALSE, // check if this is right
2605 AliAODTrack::kOrphan);
2606 aodTrack->ConvertAliPIDtoAODPID();
2608 } // end of loop on tracks
2611 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2612 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2614 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2615 p[0] = esdMuTrack->Px();
2616 p[1] = esdMuTrack->Py();
2617 p[2] = esdMuTrack->Pz();
2618 pos[0] = primary->GetX();
2619 pos[1] = primary->GetY();
2620 pos[2] = primary->GetZ();
2622 // has to be changed once the muon pid is provided by the ESD
2623 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2625 primary->AddDaughter(
2626 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2627 0, // no label provided
2632 NULL, // no covariance matrix provided
2633 (Short_t)-99, // no charge provided
2634 0, // no ITSClusterMap
2637 kTRUE, // check if this is right
2638 kTRUE, // not used for vertex fit
2639 AliAODTrack::kPrimary)
2643 // Access to the AOD container of clusters
2644 TClonesArray &clusters = *(aod->GetClusters());
2648 Int_t nClusters = esd->GetNumberOfCaloClusters();
2650 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2652 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2654 Int_t id = cluster->GetID();
2656 Float_t energy = cluster->GetClusterEnergy();
2657 cluster->GetGlobalPosition(posF);
2658 AliAODVertex *prodVertex = primary;
2659 AliAODTrack *primTrack = NULL;
2660 Char_t ttype=AliAODCluster::kUndef;
2662 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2663 else if (cluster->IsEMCAL()) {
2665 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2666 ttype = AliAODCluster::kEMCALPseudoCluster;
2668 ttype = AliAODCluster::kEMCALClusterv1;
2672 new(clusters[jClusters++]) AliAODCluster(id,
2676 NULL, // no covariance matrix provided
2677 NULL, // no pid for clusters provided
2682 } // end of loop on calo clusters
2684 delete [] usedTrack;
2688 // fill the tree for this event
2690 } // end of event loop
2692 aodTree->GetUserInfo()->Add(aod);
2697 // write the tree to the specified file
2698 aodFile = aodTree->GetCurrentFile();
2706 void AliReconstruction::WriteAlignmentData(AliESD* esd)
2708 // Write space-points which are then used in the alignment procedures
2709 // For the moment only ITS, TRD and TPC
2711 // Load TOF clusters
2713 fLoader[3]->LoadRecPoints("read");
2714 TTree* tree = fLoader[3]->TreeR();
2716 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2719 fTracker[3]->LoadClusters(tree);
2721 Int_t ntracks = esd->GetNumberOfTracks();
2722 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2724 AliESDtrack *track = esd->GetTrack(itrack);
2727 for (Int_t iDet = 3; iDet >= 0; iDet--)
2728 nsp += track->GetNcls(iDet);
2730 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2731 track->SetTrackPointArray(sp);
2733 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2734 AliTracker *tracker = fTracker[iDet];
2735 if (!tracker) continue;
2736 Int_t nspdet = track->GetNcls(iDet);
2737 if (nspdet <= 0) continue;
2738 track->GetClusters(iDet,idx);
2742 while (isp < nspdet) {
2743 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2744 const Int_t kNTPCmax = 159;
2745 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
2746 if (!isvalid) continue;
2747 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2753 fTracker[3]->UnloadClusters();
2754 fLoader[3]->UnloadRecPoints();
2758 //_____________________________________________________________________________
2759 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2761 // The method reads the raw-data error log
2762 // accumulated within the rawReader.
2763 // It extracts the raw-data errors related to
2764 // the current event and stores them into
2765 // a TClonesArray inside the esd object.
2767 if (!fRawReader) return;
2769 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2771 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2773 if (iEvent != log->GetEventNumber()) continue;
2775 esd->AddRawDataErrorLog(log);
2780 TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2782 in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2783 Int_t kBytes = (Int_t)in.tellg();
2784 printf("Size: %d \n",kBytes);
2787 char* memblock = new char [kBytes];
2788 in.seekg (0, ios::beg);
2789 in.read (memblock, kBytes);
2791 TString fData(memblock,kBytes);
2792 fn = new TNamed(fName,fData);
2793 printf("fData Size: %d \n",fData.Sizeof());
2794 printf("fName Size: %d \n",fName.Sizeof());
2795 printf("fn Size: %d \n",fn->Sizeof());
2799 AliInfo(Form("Could not Open %s\n",fPath.Data()));
2805 void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2807 // This is not really needed in AliReconstruction at the moment
2808 // but can serve as a template
2810 TList *fList = fTree->GetUserInfo();
2811 TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2812 printf("fn Size: %d \n",fn->Sizeof());
2814 TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2815 const char* cdata = fn->GetTitle();
2816 printf("fTmp Size %d\n",fTmp.Sizeof());
2818 int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2819 printf("calculated size %d\n",size);
2820 ofstream out(fName.Data(),ios::out | ios::binary);
2821 out.write(cdata,size);