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=0;
737 Bool_t dovertex=kTRUE;
738 TObject* obj = fOptions.FindObject("ITS");
740 TString optITS = obj->GetTitle();
741 if (optITS.Contains("cosmics") || optITS.Contains("COSMICS"))
744 if(dovertex) pvtx=tVertexer.FindPrimaryVertex(esd);
745 if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
748 if (pvtx->GetStatus()) {
749 // Store the improved primary vertex
750 esd->SetPrimaryVertex(pvtx);
751 // Propagate the tracks to the DCA to the improved primary vertex
752 Double_t somethingbig = 777.;
753 Double_t bz = esd->GetMagneticField();
754 Int_t nt=esd->GetNumberOfTracks();
756 AliESDtrack *t = esd->GetTrack(nt);
757 t->RelateToVertex(pvtx, bz, somethingbig);
764 vtxer.Tracks2V0vertices(esd);
767 AliCascadeVertexer cvtxer;
768 cvtxer.V0sTracks2CascadeVertices(esd);
772 if (fWriteESDfriend) {
773 new (esdf) AliESDfriend(); // Reset...
774 esd->GetESDfriend(esdf);
781 if (fCheckPointLevel > 0) WriteESD(esd, "final");
785 // delete esdf; esdf = 0;
789 tree->GetUserInfo()->Add(esd);
790 hlttree->GetUserInfo()->Add(hltesd);
794 if(fESDPar.Contains("ESD.par")){
795 AliInfo("Attaching ESD.par to Tree");
796 TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
797 tree->GetUserInfo()->Add(fn);
801 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
802 stopwatch.RealTime(),stopwatch.CpuTime()));
806 tree->SetBranchStatus("ESDfriend*",0);
811 TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
812 ESDFile2AODFile(file, aodFile);
817 // Create tags for the events in the ESD tree (the ESD tree is always present)
818 // In case of empty events the tags will contain dummy values
821 CleanUp(file, fileOld);
827 //_____________________________________________________________________________
828 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
830 // run the local reconstruction
832 TStopwatch stopwatch;
835 AliCDBManager* man = AliCDBManager::Instance();
836 Bool_t origCache = man->GetCacheFlag();
838 TString detStr = detectors;
839 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
840 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
841 AliReconstructor* reconstructor = GetReconstructor(iDet);
842 if (!reconstructor) continue;
843 if (reconstructor->HasLocalReconstruction()) continue;
845 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
846 TStopwatch stopwatchDet;
847 stopwatchDet.Start();
849 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
851 man->SetCacheFlag(kTRUE);
852 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
853 man->GetAll(calibPath); // entries are cached!
856 fRawReader->RewindEvents();
857 reconstructor->Reconstruct(fRunLoader, fRawReader);
859 reconstructor->Reconstruct(fRunLoader);
861 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
862 fgkDetectorName[iDet],
863 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
865 // unload calibration data
869 man->SetCacheFlag(origCache);
871 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
872 AliError(Form("the following detectors were not found: %s",
874 if (fStopOnError) return kFALSE;
877 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
878 stopwatch.RealTime(),stopwatch.CpuTime()));
883 //_____________________________________________________________________________
884 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
886 // run the local reconstruction
888 TStopwatch stopwatch;
891 TString detStr = detectors;
892 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
893 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
894 AliReconstructor* reconstructor = GetReconstructor(iDet);
895 if (!reconstructor) continue;
896 AliLoader* loader = fLoader[iDet];
898 // conversion of digits
899 if (fRawReader && reconstructor->HasDigitConversion()) {
900 AliInfo(Form("converting raw data digits into root objects for %s",
901 fgkDetectorName[iDet]));
902 TStopwatch stopwatchDet;
903 stopwatchDet.Start();
904 loader->LoadDigits("update");
905 loader->CleanDigits();
906 loader->MakeDigitsContainer();
907 TTree* digitsTree = loader->TreeD();
908 reconstructor->ConvertDigits(fRawReader, digitsTree);
909 loader->WriteDigits("OVERWRITE");
910 loader->UnloadDigits();
911 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
912 fgkDetectorName[iDet],
913 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
916 // local reconstruction
917 if (!reconstructor->HasLocalReconstruction()) continue;
918 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
919 TStopwatch stopwatchDet;
920 stopwatchDet.Start();
921 loader->LoadRecPoints("update");
922 loader->CleanRecPoints();
923 loader->MakeRecPointsContainer();
924 TTree* clustersTree = loader->TreeR();
925 if (fRawReader && !reconstructor->HasDigitConversion()) {
926 reconstructor->Reconstruct(fRawReader, clustersTree);
928 loader->LoadDigits("read");
929 TTree* digitsTree = loader->TreeD();
931 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
932 if (fStopOnError) return kFALSE;
934 reconstructor->Reconstruct(digitsTree, clustersTree);
936 loader->UnloadDigits();
938 loader->WriteRecPoints("OVERWRITE");
939 loader->UnloadRecPoints();
940 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
941 fgkDetectorName[iDet],
942 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
945 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
946 AliError(Form("the following detectors were not found: %s",
948 if (fStopOnError) return kFALSE;
951 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
952 stopwatch.RealTime(),stopwatch.CpuTime()));
957 //_____________________________________________________________________________
958 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
960 // run the barrel tracking
962 TStopwatch stopwatch;
965 AliESDVertex* vertex = NULL;
966 Double_t vtxPos[3] = {0, 0, 0};
967 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
969 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
970 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
971 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
975 if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
976 AliInfo("running the ITS vertex finder");
977 if (fLoader[0]) fLoader[0]->LoadRecPoints();
978 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
979 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
981 AliWarning("Vertex not found");
982 vertex = new AliESDVertex();
983 vertex->SetName("default");
986 vertex->SetName("reconstructed");
990 AliInfo("getting the primary vertex from MC");
991 vertex = new AliESDVertex(vtxPos, vtxErr);
995 vertex->GetXYZ(vtxPos);
996 vertex->GetSigmaXYZ(vtxErr);
998 AliWarning("no vertex reconstructed");
999 vertex = new AliESDVertex(vtxPos, vtxErr);
1001 esd->SetVertex(vertex);
1002 // if SPD multiplicity has been determined, it is stored in the ESD
1003 AliMultiplicity *mult = fVertexer->GetMultiplicity();
1004 if(mult)esd->SetMultiplicity(mult);
1006 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1007 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1011 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1012 stopwatch.RealTime(),stopwatch.CpuTime()));
1017 //_____________________________________________________________________________
1018 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1020 // run the HLT barrel tracking
1022 TStopwatch stopwatch;
1026 AliError("Missing runLoader!");
1030 AliInfo("running HLT tracking");
1032 // Get a pointer to the HLT reconstructor
1033 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1034 if (!reconstructor) return kFALSE;
1037 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1038 TString detName = fgkDetectorName[iDet];
1039 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1040 reconstructor->SetOption(detName.Data());
1041 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1043 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1044 if (fStopOnError) return kFALSE;
1048 Double_t vtxErr[3]={0.005,0.005,0.010};
1049 const AliESDVertex *vertex = esd->GetVertex();
1050 vertex->GetXYZ(vtxPos);
1051 tracker->SetVertex(vtxPos,vtxErr);
1053 fLoader[iDet]->LoadRecPoints("read");
1054 TTree* tree = fLoader[iDet]->TreeR();
1056 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1059 tracker->LoadClusters(tree);
1061 if (tracker->Clusters2Tracks(esd) != 0) {
1062 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1066 tracker->UnloadClusters();
1071 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1072 stopwatch.RealTime(),stopwatch.CpuTime()));
1077 //_____________________________________________________________________________
1078 Bool_t AliReconstruction::RunMuonTracking()
1080 // run the muon spectrometer tracking
1082 TStopwatch stopwatch;
1086 AliError("Missing runLoader!");
1089 Int_t iDet = 7; // for MUON
1091 AliInfo("is running...");
1093 // Get a pointer to the MUON reconstructor
1094 AliReconstructor *reconstructor = GetReconstructor(iDet);
1095 if (!reconstructor) return kFALSE;
1098 TString detName = fgkDetectorName[iDet];
1099 AliDebug(1, Form("%s tracking", detName.Data()));
1100 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1102 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1107 fLoader[iDet]->LoadTracks("update");
1108 fLoader[iDet]->CleanTracks();
1109 fLoader[iDet]->MakeTracksContainer();
1112 fLoader[iDet]->LoadRecPoints("read");
1114 if (!tracker->Clusters2Tracks(0x0)) {
1115 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1118 fLoader[iDet]->UnloadRecPoints();
1120 fLoader[iDet]->WriteTracks("OVERWRITE");
1121 fLoader[iDet]->UnloadTracks();
1126 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1127 stopwatch.RealTime(),stopwatch.CpuTime()));
1133 //_____________________________________________________________________________
1134 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1136 // run the barrel tracking
1138 TStopwatch stopwatch;
1141 AliInfo("running tracking");
1143 //Fill the ESD with the T0 info (will be used by the TOF)
1144 if (fReconstructor[11])
1145 GetReconstructor(11)->FillESD(fRunLoader, esd);
1147 // pass 1: TPC + ITS inwards
1148 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1149 if (!fTracker[iDet]) continue;
1150 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1153 fLoader[iDet]->LoadRecPoints("read");
1154 TTree* tree = fLoader[iDet]->TreeR();
1156 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1159 fTracker[iDet]->LoadClusters(tree);
1162 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1163 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1166 if (fCheckPointLevel > 1) {
1167 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1169 // preliminary PID in TPC needed by the ITS tracker
1171 GetReconstructor(1)->FillESD(fRunLoader, esd);
1172 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1173 AliESDpid::MakePID(esd);
1177 // pass 2: ALL backwards
1178 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1179 if (!fTracker[iDet]) continue;
1180 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1183 if (iDet > 1) { // all except ITS, TPC
1185 fLoader[iDet]->LoadRecPoints("read");
1186 tree = fLoader[iDet]->TreeR();
1188 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1191 fTracker[iDet]->LoadClusters(tree);
1195 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1196 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1199 if (fCheckPointLevel > 1) {
1200 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1204 if (iDet > 2) { // all except ITS, TPC, TRD
1205 fTracker[iDet]->UnloadClusters();
1206 fLoader[iDet]->UnloadRecPoints();
1208 // updated PID in TPC needed by the ITS tracker -MI
1210 GetReconstructor(1)->FillESD(fRunLoader, esd);
1211 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1212 AliESDpid::MakePID(esd);
1216 // write space-points to the ESD in case alignment data output
1218 if (fWriteAlignmentData)
1219 WriteAlignmentData(esd);
1221 // pass 3: TRD + TPC + ITS refit inwards
1222 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1223 if (!fTracker[iDet]) continue;
1224 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1227 if (fTracker[iDet]->RefitInward(esd) != 0) {
1228 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1231 if (fCheckPointLevel > 1) {
1232 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1236 fTracker[iDet]->UnloadClusters();
1237 fLoader[iDet]->UnloadRecPoints();
1240 // Propagate track to the vertex - if not done by ITS
1242 Int_t ntracks = esd->GetNumberOfTracks();
1243 for (Int_t itrack=0; itrack<ntracks; itrack++){
1244 const Double_t kRadius = 3; // beam pipe radius
1245 const Double_t kMaxStep = 5; // max step
1246 const Double_t kMaxD = 123456; // max distance to prim vertex
1247 Double_t fieldZ = AliTracker::GetBz(); //
1248 AliESDtrack * track = esd->GetTrack(itrack);
1249 if (!track) continue;
1250 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1251 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1252 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1255 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1256 stopwatch.RealTime(),stopwatch.CpuTime()));
1261 //_____________________________________________________________________________
1262 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1264 // fill the event summary data
1266 TStopwatch stopwatch;
1268 AliInfo("filling ESD");
1270 TString detStr = detectors;
1271 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1272 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1273 AliReconstructor* reconstructor = GetReconstructor(iDet);
1274 if (!reconstructor) continue;
1276 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1277 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1278 TTree* clustersTree = NULL;
1279 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1280 fLoader[iDet]->LoadRecPoints("read");
1281 clustersTree = fLoader[iDet]->TreeR();
1282 if (!clustersTree) {
1283 AliError(Form("Can't get the %s clusters tree",
1284 fgkDetectorName[iDet]));
1285 if (fStopOnError) return kFALSE;
1288 if (fRawReader && !reconstructor->HasDigitConversion()) {
1289 reconstructor->FillESD(fRawReader, clustersTree, esd);
1291 TTree* digitsTree = NULL;
1292 if (fLoader[iDet]) {
1293 fLoader[iDet]->LoadDigits("read");
1294 digitsTree = fLoader[iDet]->TreeD();
1296 AliError(Form("Can't get the %s digits tree",
1297 fgkDetectorName[iDet]));
1298 if (fStopOnError) return kFALSE;
1301 reconstructor->FillESD(digitsTree, clustersTree, esd);
1302 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1304 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1305 fLoader[iDet]->UnloadRecPoints();
1309 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1311 reconstructor->FillESD(fRunLoader, esd);
1313 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1317 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1318 AliError(Form("the following detectors were not found: %s",
1320 if (fStopOnError) return kFALSE;
1323 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1324 stopwatch.RealTime(),stopwatch.CpuTime()));
1329 //_____________________________________________________________________________
1330 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1332 // Reads the trigger decision which is
1333 // stored in Trigger.root file and fills
1334 // the corresponding esd entries
1336 AliInfo("Filling trigger information into the ESD");
1339 AliCTPRawStream input(fRawReader);
1340 if (!input.Next()) {
1341 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1344 esd->SetTriggerMask(input.GetClassMask());
1345 esd->SetTriggerCluster(input.GetClusterMask());
1348 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1350 if (!runloader->LoadTrigger()) {
1351 AliCentralTrigger *aCTP = runloader->GetTrigger();
1352 esd->SetTriggerMask(aCTP->GetClassMask());
1353 esd->SetTriggerCluster(aCTP->GetClusterMask());
1356 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1361 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1373 //_____________________________________________________________________________
1374 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1377 // Filling information from RawReader Header
1380 AliInfo("Filling information from RawReader Header");
1381 esd->SetBunchCrossNumber(0);
1382 esd->SetOrbitNumber(0);
1383 esd->SetPeriodNumber(0);
1384 esd->SetTimeStamp(0);
1385 esd->SetEventType(0);
1386 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1389 const UInt_t *id = eventHeader->GetP("Id");
1390 esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1391 esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1392 esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1394 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1395 esd->SetEventType((eventHeader->Get("Type")));
1402 //_____________________________________________________________________________
1403 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1405 // check whether detName is contained in detectors
1406 // if yes, it is removed from detectors
1408 // check if all detectors are selected
1409 if ((detectors.CompareTo("ALL") == 0) ||
1410 detectors.BeginsWith("ALL ") ||
1411 detectors.EndsWith(" ALL") ||
1412 detectors.Contains(" ALL ")) {
1417 // search for the given detector
1418 Bool_t result = kFALSE;
1419 if ((detectors.CompareTo(detName) == 0) ||
1420 detectors.BeginsWith(detName+" ") ||
1421 detectors.EndsWith(" "+detName) ||
1422 detectors.Contains(" "+detName+" ")) {
1423 detectors.ReplaceAll(detName, "");
1427 // clean up the detectors string
1428 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1429 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1430 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1435 //_____________________________________________________________________________
1436 Bool_t AliReconstruction::InitRunLoader()
1438 // get or create the run loader
1440 if (gAlice) delete gAlice;
1443 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1444 // load all base libraries to get the loader classes
1445 TString libs = gSystem->GetLibraries();
1446 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1447 TString detName = fgkDetectorName[iDet];
1448 if (detName == "HLT") continue;
1449 if (libs.Contains("lib" + detName + "base.so")) continue;
1450 gSystem->Load("lib" + detName + "base.so");
1452 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1454 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1458 fRunLoader->CdGAFile();
1459 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1460 if (fRunLoader->LoadgAlice() == 0) {
1461 gAlice = fRunLoader->GetAliRun();
1462 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1465 if (!gAlice && !fRawReader) {
1466 AliError(Form("no gAlice object found in file %s",
1467 fGAliceFileName.Data()));
1472 } else { // galice.root does not exist
1474 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1478 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1479 AliConfig::GetDefaultEventFolderName(),
1482 AliError(Form("could not create run loader in file %s",
1483 fGAliceFileName.Data()));
1487 fRunLoader->MakeTree("E");
1489 while (fRawReader->NextEvent()) {
1490 fRunLoader->SetEventNumber(iEvent);
1491 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1493 fRunLoader->MakeTree("H");
1494 fRunLoader->TreeE()->Fill();
1497 fRawReader->RewindEvents();
1498 if (fNumberOfEventsPerFile > 0)
1499 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1501 fRunLoader->SetNumberOfEventsPerFile(iEvent);
1502 fRunLoader->WriteHeader("OVERWRITE");
1503 fRunLoader->CdGAFile();
1504 fRunLoader->Write(0, TObject::kOverwrite);
1505 // AliTracker::SetFieldMap(???);
1511 //_____________________________________________________________________________
1512 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1514 // get the reconstructor object and the loader for a detector
1516 if (fReconstructor[iDet]) return fReconstructor[iDet];
1518 // load the reconstructor object
1519 TPluginManager* pluginManager = gROOT->GetPluginManager();
1520 TString detName = fgkDetectorName[iDet];
1521 TString recName = "Ali" + detName + "Reconstructor";
1522 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1524 if (detName == "HLT") {
1525 if (!gROOT->GetClass("AliLevel3")) {
1526 gSystem->Load("libAliHLTSrc.so");
1527 gSystem->Load("libAliHLTMisc.so");
1528 gSystem->Load("libAliHLTHough.so");
1529 gSystem->Load("libAliHLTComp.so");
1533 AliReconstructor* reconstructor = NULL;
1534 // first check if a plugin is defined for the reconstructor
1535 TPluginHandler* pluginHandler =
1536 pluginManager->FindHandler("AliReconstructor", detName);
1537 // if not, add a plugin for it
1538 if (!pluginHandler) {
1539 AliDebug(1, Form("defining plugin for %s", recName.Data()));
1540 TString libs = gSystem->GetLibraries();
1541 if (libs.Contains("lib" + detName + "base.so") ||
1542 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1543 pluginManager->AddHandler("AliReconstructor", detName,
1544 recName, detName + "rec", recName + "()");
1546 pluginManager->AddHandler("AliReconstructor", detName,
1547 recName, detName, recName + "()");
1549 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1551 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1552 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1554 if (reconstructor) {
1555 TObject* obj = fOptions.FindObject(detName.Data());
1556 if (obj) reconstructor->SetOption(obj->GetTitle());
1557 reconstructor->Init(fRunLoader);
1558 fReconstructor[iDet] = reconstructor;
1561 // get or create the loader
1562 if (detName != "HLT") {
1563 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1564 if (!fLoader[iDet]) {
1565 AliConfig::Instance()
1566 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1568 // first check if a plugin is defined for the loader
1570 pluginManager->FindHandler("AliLoader", detName);
1571 // if not, add a plugin for it
1572 if (!pluginHandler) {
1573 TString loaderName = "Ali" + detName + "Loader";
1574 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1575 pluginManager->AddHandler("AliLoader", detName,
1576 loaderName, detName + "base",
1577 loaderName + "(const char*, TFolder*)");
1578 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1580 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1582 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1583 fRunLoader->GetEventFolder());
1585 if (!fLoader[iDet]) { // use default loader
1586 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1588 if (!fLoader[iDet]) {
1589 AliWarning(Form("couldn't get loader for %s", detName.Data()));
1590 if (fStopOnError) return NULL;
1592 fRunLoader->AddLoader(fLoader[iDet]);
1593 fRunLoader->CdGAFile();
1594 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1595 fRunLoader->Write(0, TObject::kOverwrite);
1600 return reconstructor;
1603 //_____________________________________________________________________________
1604 Bool_t AliReconstruction::CreateVertexer()
1606 // create the vertexer
1609 AliReconstructor* itsReconstructor = GetReconstructor(0);
1610 if (itsReconstructor) {
1611 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1614 AliWarning("couldn't create a vertexer for ITS");
1615 if (fStopOnError) return kFALSE;
1621 //_____________________________________________________________________________
1622 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1624 // create the trackers
1626 TString detStr = detectors;
1627 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1628 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1629 AliReconstructor* reconstructor = GetReconstructor(iDet);
1630 if (!reconstructor) continue;
1631 TString detName = fgkDetectorName[iDet];
1632 if (detName == "HLT") {
1633 fRunHLTTracking = kTRUE;
1636 if (detName == "MUON") {
1637 fRunMuonTracking = kTRUE;
1642 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1643 if (!fTracker[iDet] && (iDet < 7)) {
1644 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1645 if (fStopOnError) return kFALSE;
1652 //_____________________________________________________________________________
1653 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1655 // delete trackers and the run loader and close and delete the file
1657 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1658 delete fReconstructor[iDet];
1659 fReconstructor[iDet] = NULL;
1660 fLoader[iDet] = NULL;
1661 delete fTracker[iDet];
1662 fTracker[iDet] = NULL;
1666 delete fDiamondProfile;
1667 fDiamondProfile = NULL;
1682 gSystem->Unlink("AliESDs.old.root");
1687 //_____________________________________________________________________________
1688 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1690 // read the ESD event from a file
1692 if (!esd) return kFALSE;
1694 sprintf(fileName, "ESD_%d.%d_%s.root",
1695 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1696 if (gSystem->AccessPathName(fileName)) return kFALSE;
1698 AliInfo(Form("reading ESD from file %s", fileName));
1699 AliDebug(1, Form("reading ESD from file %s", fileName));
1700 TFile* file = TFile::Open(fileName);
1701 if (!file || !file->IsOpen()) {
1702 AliError(Form("opening %s failed", fileName));
1709 esd = (AliESD*) file->Get("ESD");
1715 //_____________________________________________________________________________
1716 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1718 // write the ESD event to a file
1722 sprintf(fileName, "ESD_%d.%d_%s.root",
1723 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1725 AliDebug(1, Form("writing ESD to file %s", fileName));
1726 TFile* file = TFile::Open(fileName, "recreate");
1727 if (!file || !file->IsOpen()) {
1728 AliError(Form("opening %s failed", fileName));
1739 //_____________________________________________________________________________
1740 void AliReconstruction::CreateTag(TFile* file)
1743 Float_t lhcLuminosity = 0.0;
1744 TString lhcState = "test";
1745 UInt_t detectorMask = 0;
1750 Double_t fMUONMASS = 0.105658369;
1753 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1754 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1756 TLorentzVector fEPvector;
1758 Float_t fZVertexCut = 10.0;
1759 Float_t fRhoVertexCut = 2.0;
1761 Float_t fLowPtCut = 1.0;
1762 Float_t fHighPtCut = 3.0;
1763 Float_t fVeryHighPtCut = 10.0;
1766 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1768 // Creates the tags for all the events in a given ESD file
1770 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1771 Int_t nPos, nNeg, nNeutr;
1772 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1773 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1774 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1775 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1776 Float_t maxPt = .0, meanPt = .0, totalP = .0;
1778 Int_t iRunNumber = 0;
1779 TString fVertexName("default");
1781 AliRunTag *tag = new AliRunTag();
1782 AliEventTag *evTag = new AliEventTag();
1783 TTree ttag("T","A Tree with event tags");
1784 TBranch * btag = ttag.Branch("AliTAG", &tag);
1785 btag->SetCompressionLevel(9);
1787 AliInfo(Form("Creating the tags......."));
1789 if (!file || !file->IsOpen()) {
1790 AliError(Form("opening failed"));
1794 Int_t lastEvent = 0;
1795 TTree *b = (TTree*) file->Get("esdTree");
1796 AliESD *esd = new AliESD();
1797 esd->ReadFromTree(b);
1799 b->GetEntry(fFirstEvent);
1800 Int_t iInitRunNumber = esd->GetRunNumber();
1802 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1803 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1804 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1832 b->GetEntry(iEventNumber);
1833 iRunNumber = esd->GetRunNumber();
1834 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1835 const AliESDVertex * vertexIn = esd->GetVertex();
1836 if (!vertexIn) AliError("ESD has not defined vertex.");
1837 if (vertexIn) fVertexName = vertexIn->GetName();
1838 if(fVertexName != "default") fVertexflag = 1;
1839 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1840 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1841 UInt_t status = esdTrack->GetStatus();
1843 //select only tracks with ITS refit
1844 if ((status&AliESDtrack::kITSrefit)==0) continue;
1845 //select only tracks with TPC refit
1846 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1848 //select only tracks with the "combined PID"
1849 if ((status&AliESDtrack::kESDpid)==0) continue;
1851 esdTrack->GetPxPyPz(p);
1852 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1853 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1856 if(fPt > maxPt) maxPt = fPt;
1858 if(esdTrack->GetSign() > 0) {
1860 if(fPt > fLowPtCut) nCh1GeV++;
1861 if(fPt > fHighPtCut) nCh3GeV++;
1862 if(fPt > fVeryHighPtCut) nCh10GeV++;
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) nNeutr++;
1874 esdTrack->GetESDpid(prob);
1877 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1878 if(rcc == 0.0) continue;
1881 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1884 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1886 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1888 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1890 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1892 if(fPt > fLowPtCut) nEl1GeV++;
1893 if(fPt > fHighPtCut) nEl3GeV++;
1894 if(fPt > fVeryHighPtCut) nEl10GeV++;
1902 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1903 // loop over all reconstructed tracks (also first track of combination)
1904 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1905 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1906 if (muonTrack == 0x0) continue;
1908 // Coordinates at vertex
1909 fZ = muonTrack->GetZ();
1910 fY = muonTrack->GetBendingCoor();
1911 fX = muonTrack->GetNonBendingCoor();
1913 fThetaX = muonTrack->GetThetaX();
1914 fThetaY = muonTrack->GetThetaY();
1916 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1917 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1918 fPxRec = fPzRec * TMath::Tan(fThetaX);
1919 fPyRec = fPzRec * TMath::Tan(fThetaY);
1920 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1922 //ChiSquare of the track if needed
1923 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1924 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1925 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1927 // total number of muons inside a vertex cut
1928 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1930 if(fEPvector.Pt() > fLowPtCut) {
1932 if(fEPvector.Pt() > fHighPtCut) {
1934 if (fEPvector.Pt() > fVeryHighPtCut) {
1942 // Fill the event tags
1944 meanPt = meanPt/ntrack;
1946 evTag->SetEventId(iEventNumber+1);
1948 evTag->SetVertexX(vertexIn->GetXv());
1949 evTag->SetVertexY(vertexIn->GetYv());
1950 evTag->SetVertexZ(vertexIn->GetZv());
1951 evTag->SetVertexZError(vertexIn->GetZRes());
1953 evTag->SetVertexFlag(fVertexflag);
1955 evTag->SetT0VertexZ(esd->GetT0zVertex());
1957 evTag->SetTriggerMask(esd->GetTriggerMask());
1958 evTag->SetTriggerCluster(esd->GetTriggerCluster());
1960 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1961 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1962 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1963 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1964 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1965 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1968 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1969 evTag->SetNumOfPosTracks(nPos);
1970 evTag->SetNumOfNegTracks(nNeg);
1971 evTag->SetNumOfNeutrTracks(nNeutr);
1973 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1974 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1975 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1976 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1978 evTag->SetNumOfProtons(nProtons);
1979 evTag->SetNumOfKaons(nKaons);
1980 evTag->SetNumOfPions(nPions);
1981 evTag->SetNumOfMuons(nMuons);
1982 evTag->SetNumOfElectrons(nElectrons);
1983 evTag->SetNumOfPhotons(nGamas);
1984 evTag->SetNumOfPi0s(nPi0s);
1985 evTag->SetNumOfNeutrons(nNeutrons);
1986 evTag->SetNumOfKaon0s(nK0s);
1988 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1989 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1990 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1991 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1992 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1993 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1994 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1995 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1996 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1998 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1999 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2001 evTag->SetTotalMomentum(totalP);
2002 evTag->SetMeanPt(meanPt);
2003 evTag->SetMaxPt(maxPt);
2005 tag->SetLHCTag(lhcLuminosity,lhcState);
2006 tag->SetDetectorTag(detectorMask);
2008 tag->SetRunId(iInitRunNumber);
2009 tag->AddEventTag(*evTag);
2011 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2012 else lastEvent = fLastEvent;
2018 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2019 tag->GetRunId(),fFirstEvent,lastEvent );
2020 AliInfo(Form("writing tags to file %s", fileName));
2021 AliDebug(1, Form("writing tags to file %s", fileName));
2023 TFile* ftag = TFile::Open(fileName, "recreate");
2032 //_____________________________________________________________________________
2033 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
2035 // write all files from the given esd file to an aod file
2037 // create an AliAOD object
2038 AliAODEvent *aod = new AliAODEvent();
2039 aod->CreateStdContent();
2045 TTree *aodTree = new TTree("AOD", "AliAOD tree");
2046 aodTree->Branch(aod->GetList());
2049 TTree *t = (TTree*) esdFile->Get("esdTree");
2050 TBranch *b = t->GetBranch("ESD");
2052 b->SetAddress(&esd);
2054 Int_t nEvents = b->GetEntries();
2056 // set arrays and pointers
2064 // loop over events and fill them
2065 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2066 b->GetEntry(iEvent);
2068 // Multiplicity information needed by the header (to be revised!)
2069 Int_t nTracks = esd->GetNumberOfTracks();
2070 Int_t nPosTracks = 0;
2071 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2072 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2074 // Update the header
2075 AliAODHeader* header = aod->GetHeader();
2077 header->SetRunNumber (esd->GetRunNumber() );
2078 header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
2079 header->SetOrbitNumber (esd->GetOrbitNumber() );
2080 header->SetPeriodNumber (esd->GetPeriodNumber() );
2081 header->SetTriggerMask (esd->GetTriggerMask() );
2082 header->SetTriggerCluster (esd->GetTriggerCluster() );
2083 header->SetEventType (esd->GetEventType() );
2084 header->SetMagneticField (esd->GetMagneticField() );
2085 header->SetZDCN1Energy (esd->GetZDCN1Energy() );
2086 header->SetZDCP1Energy (esd->GetZDCP1Energy() );
2087 header->SetZDCN2Energy (esd->GetZDCN2Energy() );
2088 header->SetZDCP2Energy (esd->GetZDCP2Energy() );
2089 header->SetZDCEMEnergy (esd->GetZDCEMEnergy() );
2090 header->SetRefMultiplicity (nTracks);
2091 header->SetRefMultiplicityPos(nPosTracks);
2092 header->SetRefMultiplicityNeg(nTracks - nPosTracks);
2093 header->SetMuonMagFieldScale(-999.); // FIXME
2094 header->SetCentrality(-999.); // FIXME
2098 Int_t nV0s = esd->GetNumberOfV0s();
2099 Int_t nCascades = esd->GetNumberOfCascades();
2100 Int_t nKinks = esd->GetNumberOfKinks();
2101 Int_t nVertices = nV0s + nCascades + nKinks;
2103 aod->ResetStd(nTracks, nVertices);
2104 AliAODTrack *aodTrack;
2107 // Array to take into account the tracks already added to the AOD
2108 Bool_t * usedTrack = NULL;
2110 usedTrack = new Bool_t[nTracks];
2111 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2113 // Array to take into account the V0s already added to the AOD
2114 Bool_t * usedV0 = NULL;
2116 usedV0 = new Bool_t[nV0s];
2117 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2119 // Array to take into account the kinks already added to the AOD
2120 Bool_t * usedKink = NULL;
2122 usedKink = new Bool_t[nKinks];
2123 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2126 // Access to the AOD container of vertices
2127 TClonesArray &vertices = *(aod->GetVertices());
2130 // Access to the AOD container of tracks
2131 TClonesArray &tracks = *(aod->GetTracks());
2134 // Add primary vertex. The primary tracks will be defined
2135 // after the loops on the composite objects (V0, cascades, kinks)
2136 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2138 vtx->GetXYZ(pos); // position
2139 vtx->GetCovMatrix(covVtx); //covariance matrix
2141 AliAODVertex * primary = new(vertices[jVertices++])
2142 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2144 // Create vertices starting from the most complex objects
2147 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2148 AliESDcascade *cascade = esd->GetCascade(nCascade);
2150 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2151 cascade->GetPosCovXi(covVtx);
2153 // Add the cascade vertex
2154 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2156 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2158 AliAODVertex::kCascade);
2160 primary->AddDaughter(vcascade);
2162 // Add the V0 from the cascade. The ESD class have to be optimized...
2163 // Now we have to search for the corresponding Vo in the list of V0s
2164 // using the indeces of the positive and negative tracks
2166 Int_t posFromV0 = cascade->GetPindex();
2167 Int_t negFromV0 = cascade->GetNindex();
2170 AliESDv0 * v0 = 0x0;
2173 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2175 v0 = esd->GetV0(iV0);
2176 Int_t posV0 = v0->GetPindex();
2177 Int_t negV0 = v0->GetNindex();
2179 if (posV0==posFromV0 && negV0==negFromV0) {
2185 AliAODVertex * vV0FromCascade = 0x0;
2187 if (indV0>-1 && !usedV0[indV0] ) {
2189 // the V0 exists in the array of V0s and is not used
2191 usedV0[indV0] = kTRUE;
2193 v0->GetXYZ(pos[0], pos[1], pos[2]);
2194 v0->GetPosCov(covVtx);
2196 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2198 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2203 // the V0 doesn't exist in the array of V0s or was used
2204 cerr << "Error: event " << iEvent << " cascade " << nCascade
2205 << " The V0 " << indV0
2206 << " doesn't exist in the array of V0s or was used!" << endl;
2208 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2209 cascade->GetPosCov(covVtx);
2211 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2213 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2216 vcascade->AddDaughter(vV0FromCascade);
2219 // Add the positive tracks from the V0
2221 if (! usedTrack[posFromV0]) {
2223 usedTrack[posFromV0] = kTRUE;
2225 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2226 esdTrack->GetPxPyPz(p);
2227 esdTrack->GetXYZ(pos);
2228 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2229 esdTrack->GetESDpid(pid);
2231 vV0FromCascade->AddDaughter(aodTrack =
2232 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2233 esdTrack->GetLabel(),
2239 (Short_t)esdTrack->GetSign(),
2240 esdTrack->GetITSClusterMap(),
2243 kTRUE, // check if this is right
2244 kFALSE, // check if this is right
2245 AliAODTrack::kSecondary)
2247 aodTrack->ConvertAliPIDtoAODPID();
2250 cerr << "Error: event " << iEvent << " cascade " << nCascade
2251 << " track " << posFromV0 << " has already been used!" << endl;
2254 // Add the negative tracks from the V0
2256 if (!usedTrack[negFromV0]) {
2258 usedTrack[negFromV0] = kTRUE;
2260 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2261 esdTrack->GetPxPyPz(p);
2262 esdTrack->GetXYZ(pos);
2263 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2264 esdTrack->GetESDpid(pid);
2266 vV0FromCascade->AddDaughter(aodTrack =
2267 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2268 esdTrack->GetLabel(),
2274 (Short_t)esdTrack->GetSign(),
2275 esdTrack->GetITSClusterMap(),
2278 kTRUE, // check if this is right
2279 kFALSE, // check if this is right
2280 AliAODTrack::kSecondary)
2282 aodTrack->ConvertAliPIDtoAODPID();
2285 cerr << "Error: event " << iEvent << " cascade " << nCascade
2286 << " track " << negFromV0 << " has already been used!" << endl;
2289 // Add the bachelor track from the cascade
2291 Int_t bachelor = cascade->GetBindex();
2293 if(!usedTrack[bachelor]) {
2295 usedTrack[bachelor] = kTRUE;
2297 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2298 esdTrack->GetPxPyPz(p);
2299 esdTrack->GetXYZ(pos);
2300 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2301 esdTrack->GetESDpid(pid);
2303 vcascade->AddDaughter(aodTrack =
2304 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2305 esdTrack->GetLabel(),
2311 (Short_t)esdTrack->GetSign(),
2312 esdTrack->GetITSClusterMap(),
2315 kTRUE, // check if this is right
2316 kFALSE, // check if this is right
2317 AliAODTrack::kSecondary)
2319 aodTrack->ConvertAliPIDtoAODPID();
2322 cerr << "Error: event " << iEvent << " cascade " << nCascade
2323 << " track " << bachelor << " has already been used!" << endl;
2326 // Add the primary track of the cascade (if any)
2328 } // end of the loop on cascades
2332 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2334 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2336 AliESDv0 *v0 = esd->GetV0(nV0);
2338 v0->GetXYZ(pos[0], pos[1], pos[2]);
2339 v0->GetPosCov(covVtx);
2341 AliAODVertex * vV0 =
2342 new(vertices[jVertices++]) AliAODVertex(pos,
2344 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2347 primary->AddDaughter(vV0);
2349 Int_t posFromV0 = v0->GetPindex();
2350 Int_t negFromV0 = v0->GetNindex();
2352 // Add the positive tracks from the V0
2354 if (!usedTrack[posFromV0]) {
2356 usedTrack[posFromV0] = kTRUE;
2358 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2359 esdTrack->GetPxPyPz(p);
2360 esdTrack->GetXYZ(pos);
2361 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2362 esdTrack->GetESDpid(pid);
2364 vV0->AddDaughter(aodTrack =
2365 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2366 esdTrack->GetLabel(),
2372 (Short_t)esdTrack->GetSign(),
2373 esdTrack->GetITSClusterMap(),
2376 kTRUE, // check if this is right
2377 kFALSE, // check if this is right
2378 AliAODTrack::kSecondary)
2380 aodTrack->ConvertAliPIDtoAODPID();
2383 cerr << "Error: event " << iEvent << " V0 " << nV0
2384 << " track " << posFromV0 << " has already been used!" << endl;
2387 // Add the negative tracks from the V0
2389 if (!usedTrack[negFromV0]) {
2391 usedTrack[negFromV0] = kTRUE;
2393 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2394 esdTrack->GetPxPyPz(p);
2395 esdTrack->GetXYZ(pos);
2396 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2397 esdTrack->GetESDpid(pid);
2399 vV0->AddDaughter(aodTrack =
2400 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2401 esdTrack->GetLabel(),
2407 (Short_t)esdTrack->GetSign(),
2408 esdTrack->GetITSClusterMap(),
2411 kTRUE, // check if this is right
2412 kFALSE, // check if this is right
2413 AliAODTrack::kSecondary)
2415 aodTrack->ConvertAliPIDtoAODPID();
2418 cerr << "Error: event " << iEvent << " V0 " << nV0
2419 << " track " << negFromV0 << " has already been used!" << endl;
2422 } // end of the loop on V0s
2424 // Kinks: it is a big mess the access to the information in the kinks
2425 // The loop is on the tracks in order to find the mother and daugther of each kink
2428 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2431 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2433 Int_t ikink = esdTrack->GetKinkIndex(0);
2436 // Negative kink index: mother, positive: daughter
2438 // Search for the second track of the kink
2440 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2442 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2444 Int_t jkink = esdTrack1->GetKinkIndex(0);
2446 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2448 // The two tracks are from the same kink
2450 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2453 Int_t idaughter = -1;
2455 if (ikink<0 && jkink>0) {
2460 else if (ikink>0 && jkink<0) {
2466 cerr << "Error: Wrong combination of kink indexes: "
2467 << ikink << " " << jkink << endl;
2471 // Add the mother track
2473 AliAODTrack * mother = NULL;
2475 if (!usedTrack[imother]) {
2477 usedTrack[imother] = kTRUE;
2479 AliESDtrack *esdTrack = esd->GetTrack(imother);
2480 esdTrack->GetPxPyPz(p);
2481 esdTrack->GetXYZ(pos);
2482 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2483 esdTrack->GetESDpid(pid);
2486 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2487 esdTrack->GetLabel(),
2493 (Short_t)esdTrack->GetSign(),
2494 esdTrack->GetITSClusterMap(),
2497 kTRUE, // check if this is right
2498 kTRUE, // check if this is right
2499 AliAODTrack::kPrimary);
2500 primary->AddDaughter(mother);
2501 mother->ConvertAliPIDtoAODPID();
2504 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2505 << " track " << imother << " has already been used!" << endl;
2508 // Add the kink vertex
2509 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2511 AliAODVertex * vkink =
2512 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2516 AliAODVertex::kKink);
2517 // Add the daughter track
2519 AliAODTrack * daughter = NULL;
2521 if (!usedTrack[idaughter]) {
2523 usedTrack[idaughter] = kTRUE;
2525 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2526 esdTrack->GetPxPyPz(p);
2527 esdTrack->GetXYZ(pos);
2528 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2529 esdTrack->GetESDpid(pid);
2532 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2533 esdTrack->GetLabel(),
2539 (Short_t)esdTrack->GetSign(),
2540 esdTrack->GetITSClusterMap(),
2543 kTRUE, // check if this is right
2544 kTRUE, // check if this is right
2545 AliAODTrack::kPrimary);
2546 vkink->AddDaughter(daughter);
2547 daughter->ConvertAliPIDtoAODPID();
2550 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2551 << " track " << idaughter << " has already been used!" << endl;
2563 // Tracks (primary and orphan)
2565 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2568 if (usedTrack[nTrack]) continue;
2570 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2571 esdTrack->GetPxPyPz(p);
2572 esdTrack->GetXYZ(pos);
2573 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2574 esdTrack->GetESDpid(pid);
2576 Float_t impactXY, impactZ;
2578 esdTrack->GetImpactParameters(impactXY,impactZ);
2581 // track inside the beam pipe
2583 primary->AddDaughter(aodTrack =
2584 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2585 esdTrack->GetLabel(),
2591 (Short_t)esdTrack->GetSign(),
2592 esdTrack->GetITSClusterMap(),
2595 kTRUE, // check if this is right
2596 kTRUE, // check if this is right
2597 AliAODTrack::kPrimary)
2599 aodTrack->ConvertAliPIDtoAODPID();
2602 // outside the beam pipe: orphan track
2604 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2605 esdTrack->GetLabel(),
2611 (Short_t)esdTrack->GetSign(),
2612 esdTrack->GetITSClusterMap(),
2615 kFALSE, // check if this is right
2616 kFALSE, // check if this is right
2617 AliAODTrack::kOrphan);
2618 aodTrack->ConvertAliPIDtoAODPID();
2620 } // end of loop on tracks
2623 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2624 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2626 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2627 p[0] = esdMuTrack->Px();
2628 p[1] = esdMuTrack->Py();
2629 p[2] = esdMuTrack->Pz();
2630 pos[0] = primary->GetX();
2631 pos[1] = primary->GetY();
2632 pos[2] = primary->GetZ();
2634 // has to be changed once the muon pid is provided by the ESD
2635 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2637 primary->AddDaughter(
2638 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2639 0, // no label provided
2644 NULL, // no covariance matrix provided
2645 (Short_t)-99, // no charge provided
2646 0, // no ITSClusterMap
2649 kTRUE, // check if this is right
2650 kTRUE, // not used for vertex fit
2651 AliAODTrack::kPrimary)
2655 // Access to the AOD container of clusters
2656 TClonesArray &clusters = *(aod->GetClusters());
2660 Int_t nClusters = esd->GetNumberOfCaloClusters();
2662 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2664 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2666 Int_t id = cluster->GetID();
2668 Float_t energy = cluster->GetClusterEnergy();
2669 cluster->GetGlobalPosition(posF);
2670 AliAODVertex *prodVertex = primary;
2671 AliAODTrack *primTrack = NULL;
2672 Char_t ttype=AliAODCluster::kUndef;
2674 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2675 else if (cluster->IsEMCAL()) {
2677 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2678 ttype = AliAODCluster::kEMCALPseudoCluster;
2680 ttype = AliAODCluster::kEMCALClusterv1;
2684 new(clusters[jClusters++]) AliAODCluster(id,
2688 NULL, // no covariance matrix provided
2689 NULL, // no pid for clusters provided
2694 } // end of loop on calo clusters
2696 delete [] usedTrack;
2700 // fill the tree for this event
2702 } // end of event loop
2704 aodTree->GetUserInfo()->Add(aod);
2709 // write the tree to the specified file
2710 aodFile = aodTree->GetCurrentFile();
2717 void AliReconstruction::WriteAlignmentData(AliESD* esd)
2719 // Write space-points which are then used in the alignment procedures
2720 // For the moment only ITS, TRD and TPC
2722 // Load TOF clusters
2724 fLoader[3]->LoadRecPoints("read");
2725 TTree* tree = fLoader[3]->TreeR();
2727 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2730 fTracker[3]->LoadClusters(tree);
2732 Int_t ntracks = esd->GetNumberOfTracks();
2733 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2735 AliESDtrack *track = esd->GetTrack(itrack);
2738 for (Int_t iDet = 3; iDet >= 0; iDet--)
2739 nsp += track->GetNcls(iDet);
2741 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2742 track->SetTrackPointArray(sp);
2744 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2745 AliTracker *tracker = fTracker[iDet];
2746 if (!tracker) continue;
2747 Int_t nspdet = track->GetNcls(iDet);
2748 if (nspdet <= 0) continue;
2749 track->GetClusters(iDet,idx);
2753 while (isp < nspdet) {
2754 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2755 const Int_t kNTPCmax = 159;
2756 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
2757 if (!isvalid) continue;
2758 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2764 fTracker[3]->UnloadClusters();
2765 fLoader[3]->UnloadRecPoints();
2769 //_____________________________________________________________________________
2770 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2772 // The method reads the raw-data error log
2773 // accumulated within the rawReader.
2774 // It extracts the raw-data errors related to
2775 // the current event and stores them into
2776 // a TClonesArray inside the esd object.
2778 if (!fRawReader) return;
2780 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2782 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2784 if (iEvent != log->GetEventNumber()) continue;
2786 esd->AddRawDataErrorLog(log);
2791 TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2793 in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2794 Int_t kBytes = (Int_t)in.tellg();
2795 printf("Size: %d \n",kBytes);
2798 char* memblock = new char [kBytes];
2799 in.seekg (0, ios::beg);
2800 in.read (memblock, kBytes);
2802 TString fData(memblock,kBytes);
2803 fn = new TNamed(fName,fData);
2804 printf("fData Size: %d \n",fData.Sizeof());
2805 printf("fName Size: %d \n",fName.Sizeof());
2806 printf("fn Size: %d \n",fn->Sizeof());
2810 AliInfo(Form("Could not Open %s\n",fPath.Data()));
2816 void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2818 // This is not really needed in AliReconstruction at the moment
2819 // but can serve as a template
2821 TList *fList = fTree->GetUserInfo();
2822 TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2823 printf("fn Size: %d \n",fn->Sizeof());
2825 TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2826 const char* cdata = fn->GetTitle();
2827 printf("fTmp Size %d\n",fTmp.Sizeof());
2829 int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2830 printf("calculated size %d\n",size);
2831 ofstream out(fName.Data(),ios::out | ios::binary);
2832 out.write(cdata,size);