git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Default run number -1 (Cvetan)
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
CommitLineData
596a855f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// class for running the reconstruction //
21// //
22// Clusters and tracks are created for all detectors and all events by //
23// typing: //
24// //
25// AliReconstruction rec; //
26// rec.Run(); //
27// //
28// The Run method returns kTRUE in case of successful execution. //
c71de921 29// //
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 //
32// //
33// rec.SetInput("..."); //
34// //
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 //
40// //
b26c3770 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 //
44// //
45// rec.SetEventRange(..., ...); //
46// //
47// The index -1 (default) can be used for the last event to indicate no //
48// upper limit of the event range. //
49// //
596a855f 50// The name of the galice file can be changed from the default //
e583c30d 51// "galice.root" by passing it as argument to the AliReconstruction //
52// constructor or by //
596a855f 53// //
54// rec.SetGAliceFile("..."); //
55// //
59697224 56// The local reconstruction can be switched on or off for individual //
57// detectors by //
596a855f 58// //
59697224 59// rec.SetRunLocalReconstruction("..."); //
596a855f 60// //
61// The argument is a (case sensitive) string with the names of the //
62// detectors separated by a space. The special string "ALL" selects all //
63// available detectors. This is the default. //
64// //
c71de921 65// The reconstruction of the primary vertex position can be switched off by //
66// //
67// rec.SetRunVertexFinder(kFALSE); //
68// //
b8cd5251 69// The tracking and the creation of ESD tracks can be switched on for //
70// selected detectors by //
596a855f 71// //
b8cd5251 72// rec.SetRunTracking("..."); //
596a855f 73// //
c84a5e9e 74// Uniform/nonuniform field tracking switches (default: uniform field) //
75// //
76// rec.SetUniformFieldTracking(); ( rec.SetNonuniformFieldTracking(); ) //
77// //
596a855f 78// The filling of additional ESD information can be steered by //
79// //
80// rec.SetFillESD("..."); //
81// //
b8cd5251 82// Again, for both methods the string specifies the list of detectors. //
83// The default is "ALL". //
84// //
85// The call of the shortcut method //
86// //
87// rec.SetRunReconstruction("..."); //
88// //
89// is equivalent to calling SetRunLocalReconstruction, SetRunTracking and //
90// SetFillESD with the same detector selecting string as argument. //
596a855f 91// //
c71de921 92// The reconstruction requires digits or raw data as input. For the creation //
93// of digits and raw data have a look at the class AliSimulation. //
596a855f 94// //
24f7a148 95// For debug purposes the method SetCheckPointLevel can be used. If the //
96// argument is greater than 0, files with ESD events will be written after //
97// selected steps of the reconstruction for each event: //
98// level 1: after tracking and after filling of ESD (final) //
99// level 2: in addition after each tracking step //
100// level 3: in addition after the filling of ESD for each detector //
101// If a final check point file exists for an event, this event will be //
102// skipped in the reconstruction. The tracking and the filling of ESD for //
103// a detector will be skipped as well, if the corresponding check point //
104// file exists. The ESD event will then be loaded from the file instead. //
105// //
596a855f 106///////////////////////////////////////////////////////////////////////////////
107
024a7e64 108#include <TArrayF.h>
109#include <TFile.h>
110#include <TSystem.h>
111#include <TROOT.h>
112#include <TPluginManager.h>
fd46e2d2 113#include <TStopwatch.h>
3103d196 114#include <TGeoManager.h>
2bdb9d38 115#include <TLorentzVector.h>
596a855f 116
117#include "AliReconstruction.h"
b8cd5251 118#include "AliReconstructor.h"
815c2b38 119#include "AliLog.h"
596a855f 120#include "AliRunLoader.h"
121#include "AliRun.h"
b649205a 122#include "AliRawReaderFile.h"
123#include "AliRawReaderDate.h"
124#include "AliRawReaderRoot.h"
596a855f 125#include "AliESD.h"
2257f27e 126#include "AliESDVertex.h"
c84a5e9e 127#include "AliTracker.h"
2257f27e 128#include "AliVertexer.h"
596a855f 129#include "AliHeader.h"
130#include "AliGenEventHeader.h"
b26c3770 131#include "AliPID.h"
596a855f 132#include "AliESDpid.h"
ff8bb5ae 133#include "AliESDtrack.h"
f3a97c86 134
135#include "AliRunTag.h"
089bf903 136//#include "AliLHCTag.h"
f3a97c86 137#include "AliDetectorTag.h"
138#include "AliEventTag.h"
139
98937d93 140#include "AliTrackPointArray.h"
b0314964 141#include "AliCDBManager.h"
f3a97c86 142
596a855f 143ClassImp(AliReconstruction)
144
145
146//_____________________________________________________________________________
482070f2 147const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
c757bafd 148
149//_____________________________________________________________________________
024cf675 150AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 151 const char* name, const char* title) :
152 TNamed(name, title),
153
59697224 154 fRunLocalReconstruction("ALL"),
c84a5e9e 155 fUniformField(kTRUE),
2257f27e 156 fRunVertexFinder(kTRUE),
1f46a9ae 157 fRunHLTTracking(kFALSE),
b8cd5251 158 fRunTracking("ALL"),
e583c30d 159 fFillESD("ALL"),
160 fGAliceFileName(gAliceFilename),
b649205a 161 fInput(""),
b26c3770 162 fFirstEvent(0),
163 fLastEvent(-1),
e583c30d 164 fStopOnError(kFALSE),
24f7a148 165 fCheckPointLevel(0),
b8cd5251 166 fOptions(),
e583c30d 167
168 fRunLoader(NULL),
b649205a 169 fRawReader(NULL),
b8cd5251 170
98937d93 171 fVertexer(NULL),
172
024cf675 173 fWriteAlignmentData(kFALSE),
174 fCDBUri(cdbUri)
596a855f 175{
176// create reconstruction object with default parameters
b8cd5251 177
178 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
179 fReconstructor[iDet] = NULL;
180 fLoader[iDet] = NULL;
181 fTracker[iDet] = NULL;
182 }
e47c4c2e 183 AliPID pid;
3103d196 184 // Import TGeo geometry
af83cb85 185 TString geom(gSystem->DirName(gAliceFilename));
186 geom += "/geometry.root";
187 TGeoManager::Import(geom.Data());
596a855f 188}
189
190//_____________________________________________________________________________
191AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 192 TNamed(rec),
193
59697224 194 fRunLocalReconstruction(rec.fRunLocalReconstruction),
c84a5e9e 195 fUniformField(rec.fUniformField),
2257f27e 196 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 197 fRunHLTTracking(rec.fRunHLTTracking),
e583c30d 198 fRunTracking(rec.fRunTracking),
199 fFillESD(rec.fFillESD),
200 fGAliceFileName(rec.fGAliceFileName),
b649205a 201 fInput(rec.fInput),
b26c3770 202 fFirstEvent(rec.fFirstEvent),
203 fLastEvent(rec.fLastEvent),
e583c30d 204 fStopOnError(rec.fStopOnError),
24f7a148 205 fCheckPointLevel(0),
b8cd5251 206 fOptions(),
e583c30d 207
208 fRunLoader(NULL),
b649205a 209 fRawReader(NULL),
b8cd5251 210
98937d93 211 fVertexer(NULL),
212
024cf675 213 fWriteAlignmentData(rec.fWriteAlignmentData),
214 fCDBUri(rec.fCDBUri)
596a855f 215{
216// copy constructor
217
efd2085e 218 for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
219 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
220 }
b8cd5251 221 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
222 fReconstructor[iDet] = NULL;
223 fLoader[iDet] = NULL;
224 fTracker[iDet] = NULL;
225 }
596a855f 226}
227
228//_____________________________________________________________________________
229AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
230{
231// assignment operator
232
233 this->~AliReconstruction();
234 new(this) AliReconstruction(rec);
235 return *this;
236}
237
238//_____________________________________________________________________________
239AliReconstruction::~AliReconstruction()
240{
241// clean up
242
e583c30d 243 CleanUp();
efd2085e 244 fOptions.Delete();
596a855f 245}
246
024cf675 247//_____________________________________________________________________________
248void AliReconstruction::InitCDBStorage()
249{
250// activate a default CDB storage
251// First check if we have any CDB storage set, because it is used
252// to retrieve the calibration and alignment constants
253
254 AliCDBManager* man = AliCDBManager::Instance();
255 if (!man->IsDefaultStorageSet())
256 {
257 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
258 AliWarning("Default CDB storage not yet set");
259 AliWarning(Form("Using default storage declared in AliSimulation: %s",fCDBUri.Data()));
260 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
261 SetDefaultStorage(fCDBUri);
262 }
263
264}
265
266//_____________________________________________________________________________
267void AliReconstruction::SetDefaultStorage(const char* uri) {
268// activate a default CDB storage
269
270 AliCDBManager::Instance()->SetDefaultStorage(uri);
271
272}
273
274//_____________________________________________________________________________
275void AliReconstruction::SetSpecificStorage(const char* detName, const char* uri) {
276// activate a detector-specific CDB storage
277
278 AliCDBManager::Instance()->SetSpecificStorage(detName, uri);
279
280}
281
596a855f 282
283//_____________________________________________________________________________
284void AliReconstruction::SetGAliceFile(const char* fileName)
285{
286// set the name of the galice file
287
288 fGAliceFileName = fileName;
289}
290
efd2085e 291//_____________________________________________________________________________
292void AliReconstruction::SetOption(const char* detector, const char* option)
293{
294// set options for the reconstruction of a detector
295
296 TObject* obj = fOptions.FindObject(detector);
297 if (obj) fOptions.Remove(obj);
298 fOptions.Add(new TNamed(detector, option));
299}
300
596a855f 301
302//_____________________________________________________________________________
b26c3770 303Bool_t AliReconstruction::Run(const char* input,
304 Int_t firstEvent, Int_t lastEvent)
596a855f 305{
306// run the reconstruction
307
024cf675 308 InitCDBStorage();
b0314964 309
b649205a 310 // set the input
311 if (!input) input = fInput.Data();
312 TString fileName(input);
313 if (fileName.EndsWith("/")) {
314 fRawReader = new AliRawReaderFile(fileName);
315 } else if (fileName.EndsWith(".root")) {
316 fRawReader = new AliRawReaderRoot(fileName);
317 } else if (!fileName.IsNull()) {
318 fRawReader = new AliRawReaderDate(fileName);
319 fRawReader->SelectEvents(7);
320 }
321
f08fc9f5 322 // get the run loader
323 if (!InitRunLoader()) return kFALSE;
596a855f 324
325 // local reconstruction
59697224 326 if (!fRunLocalReconstruction.IsNull()) {
327 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 328 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 329 }
330 }
b26c3770 331// if (!fRunVertexFinder && fRunTracking.IsNull() &&
332// fFillESD.IsNull()) return kTRUE;
2257f27e 333
334 // get vertexer
335 if (fRunVertexFinder && !CreateVertexer()) {
336 if (fStopOnError) {
337 CleanUp();
338 return kFALSE;
339 }
340 }
596a855f 341
f08fc9f5 342 // get trackers
b8cd5251 343 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 344 if (fStopOnError) {
345 CleanUp();
346 return kFALSE;
347 }
596a855f 348 }
24f7a148 349
9db3a215 350
351 TStopwatch stopwatch;
352 stopwatch.Start();
353
b26c3770 354 // get the possibly already existing ESD file and tree
1f46a9ae 355 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
b26c3770 356 TFile* fileOld = NULL;
1f46a9ae 357 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 358 if (!gSystem->AccessPathName("AliESDs.root")){
359 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
360 fileOld = TFile::Open("AliESDs.old.root");
361 if (fileOld && fileOld->IsOpen()) {
362 treeOld = (TTree*) fileOld->Get("esdTree");
363 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
1f46a9ae 364 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
365 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
b26c3770 366 }
367 }
368
36711aa4 369 // create the ESD output file and tree
596a855f 370 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
371 if (!file->IsOpen()) {
815c2b38 372 AliError("opening AliESDs.root failed");
b26c3770 373 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 374 }
36711aa4 375 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
376 tree->Branch("ESD", "AliESD", &esd);
1f46a9ae 377 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
378 hlttree->Branch("ESD", "AliESD", &hltesd);
379 delete esd; delete hltesd;
380 esd = NULL; hltesd = NULL;
36711aa4 381 gROOT->cd();
596a855f 382
383 // loop over events
b649205a 384 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 385
596a855f 386 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 387 if (fRawReader) fRawReader->NextEvent();
388 if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
389 // copy old ESD to the new one
390 if (treeOld) {
391 treeOld->SetBranchAddress("ESD", &esd);
392 treeOld->GetEntry(iEvent);
393 }
394 tree->Fill();
1f46a9ae 395 if (hlttreeOld) {
396 hlttreeOld->SetBranchAddress("ESD", &hltesd);
397 hlttreeOld->GetEntry(iEvent);
398 }
399 hlttree->Fill();
b26c3770 400 continue;
401 }
402
815c2b38 403 AliInfo(Form("processing event %d", iEvent));
596a855f 404 fRunLoader->GetEvent(iEvent);
24f7a148 405
406 char fileName[256];
407 sprintf(fileName, "ESD_%d.%d_final.root",
f08fc9f5 408 fRunLoader->GetHeader()->GetRun(),
409 fRunLoader->GetHeader()->GetEventNrInRun());
24f7a148 410 if (!gSystem->AccessPathName(fileName)) continue;
411
b26c3770 412 // local reconstruction
413 if (!fRunLocalReconstruction.IsNull()) {
414 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
415 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
416 }
417 }
418
1f46a9ae 419 esd = new AliESD; hltesd = new AliESD;
f08fc9f5 420 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 421 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
f08fc9f5 422 esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
1f46a9ae 423 hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
f08fc9f5 424 if (gAlice) {
7ebb06c3 425 esd->SetMagneticField(AliTracker::GetBz());
426 hltesd->SetMagneticField(AliTracker::GetBz());
f08fc9f5 427 } else {
428 // ???
429 }
596a855f 430
2257f27e 431 // vertex finder
432 if (fRunVertexFinder) {
433 if (!ReadESD(esd, "vertex")) {
434 if (!RunVertexFinder(esd)) {
b26c3770 435 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 436 }
437 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
438 }
439 }
440
1f46a9ae 441 // HLT tracking
442 if (!fRunTracking.IsNull()) {
443 if (fRunHLTTracking) {
444 hltesd->SetVertex(esd->GetVertex());
445 if (!RunHLTTracking(hltesd)) {
446 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
447 }
448 }
449 }
450
596a855f 451 // barrel tracking
b8cd5251 452 if (!fRunTracking.IsNull()) {
24f7a148 453 if (!ReadESD(esd, "tracking")) {
454 if (!RunTracking(esd)) {
b26c3770 455 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 456 }
457 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
596a855f 458 }
459 }
460
461 // fill ESD
462 if (!fFillESD.IsNull()) {
463 if (!FillESD(esd, fFillESD)) {
b26c3770 464 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 465 }
466 }
467
468 // combined PID
469 AliESDpid::MakePID(esd);
24f7a148 470 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 471
472 // write ESD
36711aa4 473 tree->Fill();
1f46a9ae 474 // write HLT ESD
475 hlttree->Fill();
24f7a148 476
f3a97c86 477 if (fCheckPointLevel > 0) WriteESD(esd, "final");
478
1f46a9ae 479 delete esd; delete hltesd;
480 esd = NULL; hltesd = NULL;
596a855f 481 }
482
9db3a215 483 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
484 stopwatch.RealTime(),stopwatch.CpuTime()));
485
36711aa4 486 file->cd();
487 tree->Write();
1f46a9ae 488 hlttree->Write();
f3a97c86 489
490 // Create tags for the events in the ESD tree (the ESD tree is always present)
491 // In case of empty events the tags will contain dummy values
492 CreateTag(file);
b26c3770 493 CleanUp(file, fileOld);
596a855f 494
495 return kTRUE;
496}
497
498
499//_____________________________________________________________________________
59697224 500Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 501{
59697224 502// run the local reconstruction
596a855f 503
030b532d 504 TStopwatch stopwatch;
505 stopwatch.Start();
506
596a855f 507 TString detStr = detectors;
b8cd5251 508 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
509 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
510 AliReconstructor* reconstructor = GetReconstructor(iDet);
511 if (!reconstructor) continue;
b26c3770 512 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 513
514 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
515 TStopwatch stopwatchDet;
516 stopwatchDet.Start();
517 if (fRawReader) {
518 fRawReader->RewindEvents();
519 reconstructor->Reconstruct(fRunLoader, fRawReader);
520 } else {
521 reconstructor->Reconstruct(fRunLoader);
596a855f 522 }
5f8272e1 523 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
524 fgkDetectorName[iDet],
525 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
596a855f 526 }
527
528 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 529 AliError(Form("the following detectors were not found: %s",
530 detStr.Data()));
596a855f 531 if (fStopOnError) return kFALSE;
532 }
533
5f8272e1 534 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
535 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 536
596a855f 537 return kTRUE;
538}
539
540//_____________________________________________________________________________
b26c3770 541Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
542{
543// run the local reconstruction
544
545 TStopwatch stopwatch;
546 stopwatch.Start();
547
548 TString detStr = detectors;
549 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
550 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
551 AliReconstructor* reconstructor = GetReconstructor(iDet);
552 if (!reconstructor) continue;
553 AliLoader* loader = fLoader[iDet];
554
555 // conversion of digits
556 if (fRawReader && reconstructor->HasDigitConversion()) {
557 AliInfo(Form("converting raw data digits into root objects for %s",
558 fgkDetectorName[iDet]));
559 TStopwatch stopwatchDet;
560 stopwatchDet.Start();
561 loader->LoadDigits("update");
562 loader->CleanDigits();
563 loader->MakeDigitsContainer();
564 TTree* digitsTree = loader->TreeD();
565 reconstructor->ConvertDigits(fRawReader, digitsTree);
566 loader->WriteDigits("OVERWRITE");
567 loader->UnloadDigits();
5f8272e1 568 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
569 fgkDetectorName[iDet],
570 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 571 }
572
573 // local reconstruction
574 if (!reconstructor->HasLocalReconstruction()) continue;
575 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
576 TStopwatch stopwatchDet;
577 stopwatchDet.Start();
578 loader->LoadRecPoints("update");
579 loader->CleanRecPoints();
580 loader->MakeRecPointsContainer();
581 TTree* clustersTree = loader->TreeR();
582 if (fRawReader && !reconstructor->HasDigitConversion()) {
583 reconstructor->Reconstruct(fRawReader, clustersTree);
584 } else {
585 loader->LoadDigits("read");
586 TTree* digitsTree = loader->TreeD();
587 if (!digitsTree) {
588 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
589 if (fStopOnError) return kFALSE;
590 } else {
591 reconstructor->Reconstruct(digitsTree, clustersTree);
592 }
593 loader->UnloadDigits();
594 }
595 loader->WriteRecPoints("OVERWRITE");
596 loader->UnloadRecPoints();
5f8272e1 597 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
598 fgkDetectorName[iDet],
599 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 600 }
601
602 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
603 AliError(Form("the following detectors were not found: %s",
604 detStr.Data()));
605 if (fStopOnError) return kFALSE;
606 }
5f8272e1 607
608 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
609 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 610
611 return kTRUE;
612}
613
614//_____________________________________________________________________________
2257f27e 615Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 616{
617// run the barrel tracking
618
030b532d 619 TStopwatch stopwatch;
620 stopwatch.Start();
621
2257f27e 622 AliESDVertex* vertex = NULL;
623 Double_t vtxPos[3] = {0, 0, 0};
624 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
625 TArrayF mcVertex(3);
a6b0b91b 626 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
627 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
628 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
629 }
2257f27e 630
b8cd5251 631 if (fVertexer) {
815c2b38 632 AliInfo("running the ITS vertex finder");
b26c3770 633 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 634 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 635 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 636 if(!vertex){
815c2b38 637 AliWarning("Vertex not found");
c710f220 638 vertex = new AliESDVertex();
d1a50cb5 639 vertex->SetName("default");
2257f27e 640 }
641 else {
642 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 643 vertex->SetName("reconstructed");
2257f27e 644 }
645
646 } else {
815c2b38 647 AliInfo("getting the primary vertex from MC");
2257f27e 648 vertex = new AliESDVertex(vtxPos, vtxErr);
649 }
650
651 if (vertex) {
652 vertex->GetXYZ(vtxPos);
653 vertex->GetSigmaXYZ(vtxErr);
654 } else {
815c2b38 655 AliWarning("no vertex reconstructed");
2257f27e 656 vertex = new AliESDVertex(vtxPos, vtxErr);
657 }
658 esd->SetVertex(vertex);
b8cd5251 659 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
660 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
661 }
2257f27e 662 delete vertex;
663
5f8272e1 664 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
665 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 666
667 return kTRUE;
668}
669
670//_____________________________________________________________________________
1f46a9ae 671Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
672{
673// run the HLT barrel tracking
674
675 TStopwatch stopwatch;
676 stopwatch.Start();
677
678 if (!fRunLoader) {
679 AliError("Missing runLoader!");
680 return kFALSE;
681 }
682
683 AliInfo("running HLT tracking");
684
685 // Get a pointer to the HLT reconstructor
686 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
687 if (!reconstructor) return kFALSE;
688
689 // TPC + ITS
690 for (Int_t iDet = 1; iDet >= 0; iDet--) {
691 TString detName = fgkDetectorName[iDet];
692 AliDebug(1, Form("%s HLT tracking", detName.Data()));
693 reconstructor->SetOption(detName.Data());
694 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
695 if (!tracker) {
696 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
697 if (fStopOnError) return kFALSE;
9dcc06e1 698 continue;
1f46a9ae 699 }
700 Double_t vtxPos[3];
701 Double_t vtxErr[3]={0.005,0.005,0.010};
702 const AliESDVertex *vertex = esd->GetVertex();
703 vertex->GetXYZ(vtxPos);
704 tracker->SetVertex(vtxPos,vtxErr);
705 if(iDet != 1) {
706 fLoader[iDet]->LoadRecPoints("read");
707 TTree* tree = fLoader[iDet]->TreeR();
708 if (!tree) {
709 AliError(Form("Can't get the %s cluster tree", detName.Data()));
710 return kFALSE;
711 }
712 tracker->LoadClusters(tree);
713 }
714 if (tracker->Clusters2Tracks(esd) != 0) {
715 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
716 return kFALSE;
717 }
718 if(iDet != 1) {
719 tracker->UnloadClusters();
720 }
721 delete tracker;
722 }
723
5f8272e1 724 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
725 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 726
727 return kTRUE;
728}
729
730//_____________________________________________________________________________
2257f27e 731Bool_t AliReconstruction::RunTracking(AliESD*& esd)
732{
733// run the barrel tracking
734
735 TStopwatch stopwatch;
736 stopwatch.Start();
24f7a148 737
815c2b38 738 AliInfo("running tracking");
596a855f 739
b8cd5251 740 // pass 1: TPC + ITS inwards
741 for (Int_t iDet = 1; iDet >= 0; iDet--) {
742 if (!fTracker[iDet]) continue;
743 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 744
b8cd5251 745 // load clusters
746 fLoader[iDet]->LoadRecPoints("read");
747 TTree* tree = fLoader[iDet]->TreeR();
748 if (!tree) {
749 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 750 return kFALSE;
751 }
b8cd5251 752 fTracker[iDet]->LoadClusters(tree);
753
754 // run tracking
755 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
756 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 757 return kFALSE;
758 }
b8cd5251 759 if (fCheckPointLevel > 1) {
760 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
761 }
878e1fe1 762 // preliminary PID in TPC needed by the ITS tracker
763 if (iDet == 1) {
764 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 765 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 766 AliESDpid::MakePID(esd);
767 }
b8cd5251 768 }
596a855f 769
b8cd5251 770 // pass 2: ALL backwards
771 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
772 if (!fTracker[iDet]) continue;
773 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
774
775 // load clusters
776 if (iDet > 1) { // all except ITS, TPC
777 TTree* tree = NULL;
7b61cd9c 778 fLoader[iDet]->LoadRecPoints("read");
779 tree = fLoader[iDet]->TreeR();
b8cd5251 780 if (!tree) {
781 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 782 return kFALSE;
783 }
b8cd5251 784 fTracker[iDet]->LoadClusters(tree);
785 }
24f7a148 786
b8cd5251 787 // run tracking
788 if (fTracker[iDet]->PropagateBack(esd) != 0) {
789 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
790 return kFALSE;
791 }
792 if (fCheckPointLevel > 1) {
793 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
794 }
24f7a148 795
b8cd5251 796 // unload clusters
797 if (iDet > 2) { // all except ITS, TPC, TRD
798 fTracker[iDet]->UnloadClusters();
7b61cd9c 799 fLoader[iDet]->UnloadRecPoints();
b8cd5251 800 }
8f37df88 801 // updated PID in TPC needed by the ITS tracker -MI
802 if (iDet == 1) {
803 GetReconstructor(1)->FillESD(fRunLoader, esd);
804 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
805 AliESDpid::MakePID(esd);
806 }
b8cd5251 807 }
596a855f 808
98937d93 809 // write space-points to the ESD in case alignment data output
810 // is switched on
811 if (fWriteAlignmentData)
812 WriteAlignmentData(esd);
813
b8cd5251 814 // pass 3: TRD + TPC + ITS refit inwards
815 for (Int_t iDet = 2; iDet >= 0; iDet--) {
816 if (!fTracker[iDet]) continue;
817 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 818
b8cd5251 819 // run tracking
820 if (fTracker[iDet]->RefitInward(esd) != 0) {
821 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
822 return kFALSE;
823 }
824 if (fCheckPointLevel > 1) {
825 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
826 }
596a855f 827
b8cd5251 828 // unload clusters
829 fTracker[iDet]->UnloadClusters();
830 fLoader[iDet]->UnloadRecPoints();
831 }
ff8bb5ae 832 //
833 // Propagate track to the vertex - if not done by ITS
834 //
835 Int_t ntracks = esd->GetNumberOfTracks();
836 for (Int_t itrack=0; itrack<ntracks; itrack++){
837 const Double_t kRadius = 3; // beam pipe radius
838 const Double_t kMaxStep = 5; // max step
839 const Double_t kMaxD = 123456; // max distance to prim vertex
840 Double_t fieldZ = AliTracker::GetBz(); //
841 AliESDtrack * track = esd->GetTrack(itrack);
842 if (!track) continue;
843 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
844 track->PropagateTo(kRadius, track->GetMass(),kMaxStep,kTRUE);
845 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
846 }
847
5f8272e1 848 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
849 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 850
596a855f 851 return kTRUE;
852}
853
854//_____________________________________________________________________________
24f7a148 855Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 856{
857// fill the event summary data
858
030b532d 859 TStopwatch stopwatch;
860 stopwatch.Start();
815c2b38 861 AliInfo("filling ESD");
030b532d 862
596a855f 863 TString detStr = detectors;
b8cd5251 864 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
865 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
866 AliReconstructor* reconstructor = GetReconstructor(iDet);
867 if (!reconstructor) continue;
868
869 if (!ReadESD(esd, fgkDetectorName[iDet])) {
870 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 871 TTree* clustersTree = NULL;
872 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
873 fLoader[iDet]->LoadRecPoints("read");
874 clustersTree = fLoader[iDet]->TreeR();
875 if (!clustersTree) {
876 AliError(Form("Can't get the %s clusters tree",
877 fgkDetectorName[iDet]));
878 if (fStopOnError) return kFALSE;
879 }
880 }
881 if (fRawReader && !reconstructor->HasDigitConversion()) {
882 reconstructor->FillESD(fRawReader, clustersTree, esd);
883 } else {
884 TTree* digitsTree = NULL;
885 if (fLoader[iDet]) {
886 fLoader[iDet]->LoadDigits("read");
887 digitsTree = fLoader[iDet]->TreeD();
888 if (!digitsTree) {
889 AliError(Form("Can't get the %s digits tree",
890 fgkDetectorName[iDet]));
891 if (fStopOnError) return kFALSE;
892 }
893 }
894 reconstructor->FillESD(digitsTree, clustersTree, esd);
895 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
896 }
897 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
898 fLoader[iDet]->UnloadRecPoints();
899 }
900
b8cd5251 901 if (fRawReader) {
902 reconstructor->FillESD(fRunLoader, fRawReader, esd);
903 } else {
904 reconstructor->FillESD(fRunLoader, esd);
24f7a148 905 }
b8cd5251 906 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 907 }
908 }
909
910 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 911 AliError(Form("the following detectors were not found: %s",
912 detStr.Data()));
596a855f 913 if (fStopOnError) return kFALSE;
914 }
915
5f8272e1 916 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
917 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 918
596a855f 919 return kTRUE;
920}
921
922
923//_____________________________________________________________________________
924Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
925{
926// check whether detName is contained in detectors
927// if yes, it is removed from detectors
928
929 // check if all detectors are selected
930 if ((detectors.CompareTo("ALL") == 0) ||
931 detectors.BeginsWith("ALL ") ||
932 detectors.EndsWith(" ALL") ||
933 detectors.Contains(" ALL ")) {
934 detectors = "ALL";
935 return kTRUE;
936 }
937
938 // search for the given detector
939 Bool_t result = kFALSE;
940 if ((detectors.CompareTo(detName) == 0) ||
941 detectors.BeginsWith(detName+" ") ||
942 detectors.EndsWith(" "+detName) ||
943 detectors.Contains(" "+detName+" ")) {
944 detectors.ReplaceAll(detName, "");
945 result = kTRUE;
946 }
947
948 // clean up the detectors string
949 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
950 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
951 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
952
953 return result;
954}
e583c30d 955
956//_____________________________________________________________________________
f08fc9f5 957Bool_t AliReconstruction::InitRunLoader()
958{
959// get or create the run loader
960
961 if (gAlice) delete gAlice;
962 gAlice = NULL;
963
b26c3770 964 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
965 // load all base libraries to get the loader classes
966 TString libs = gSystem->GetLibraries();
967 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
968 TString detName = fgkDetectorName[iDet];
969 if (detName == "HLT") continue;
970 if (libs.Contains("lib" + detName + "base.so")) continue;
971 gSystem->Load("lib" + detName + "base.so");
972 }
f08fc9f5 973 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
974 if (!fRunLoader) {
975 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
976 CleanUp();
977 return kFALSE;
978 }
b26c3770 979 fRunLoader->CdGAFile();
980 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
981 if (fRunLoader->LoadgAlice() == 0) {
982 gAlice = fRunLoader->GetAliRun();
c84a5e9e 983 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 984 }
f08fc9f5 985 }
986 if (!gAlice && !fRawReader) {
987 AliError(Form("no gAlice object found in file %s",
988 fGAliceFileName.Data()));
989 CleanUp();
990 return kFALSE;
991 }
992
993 } else { // galice.root does not exist
994 if (!fRawReader) {
995 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
996 CleanUp();
997 return kFALSE;
998 }
999 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1000 AliConfig::GetDefaultEventFolderName(),
1001 "recreate");
1002 if (!fRunLoader) {
1003 AliError(Form("could not create run loader in file %s",
1004 fGAliceFileName.Data()));
1005 CleanUp();
1006 return kFALSE;
1007 }
1008 fRunLoader->MakeTree("E");
1009 Int_t iEvent = 0;
1010 while (fRawReader->NextEvent()) {
1011 fRunLoader->SetEventNumber(iEvent);
1012 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1013 iEvent, iEvent);
1014 fRunLoader->MakeTree("H");
1015 fRunLoader->TreeE()->Fill();
1016 iEvent++;
1017 }
1018 fRawReader->RewindEvents();
1019 fRunLoader->WriteHeader("OVERWRITE");
1020 fRunLoader->CdGAFile();
1021 fRunLoader->Write(0, TObject::kOverwrite);
1022// AliTracker::SetFieldMap(???);
1023 }
1024
1025 return kTRUE;
1026}
1027
1028//_____________________________________________________________________________
b8cd5251 1029AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1030{
f08fc9f5 1031// get the reconstructor object and the loader for a detector
c757bafd 1032
b8cd5251 1033 if (fReconstructor[iDet]) return fReconstructor[iDet];
1034
1035 // load the reconstructor object
1036 TPluginManager* pluginManager = gROOT->GetPluginManager();
1037 TString detName = fgkDetectorName[iDet];
1038 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1039 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1040
1041 if (detName == "HLT") {
1042 if (!gROOT->GetClass("AliLevel3")) {
1043 gSystem->Load("libAliL3Src.so");
1044 gSystem->Load("libAliL3Misc.so");
1045 gSystem->Load("libAliL3Hough.so");
1046 gSystem->Load("libAliL3Comp.so");
1047 }
1048 }
1049
1050 AliReconstructor* reconstructor = NULL;
1051 // first check if a plugin is defined for the reconstructor
1052 TPluginHandler* pluginHandler =
1053 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1054 // if not, add a plugin for it
1055 if (!pluginHandler) {
b8cd5251 1056 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1057 TString libs = gSystem->GetLibraries();
1058 if (libs.Contains("lib" + detName + "base.so") ||
1059 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1060 pluginManager->AddHandler("AliReconstructor", detName,
1061 recName, detName + "rec", recName + "()");
1062 } else {
1063 pluginManager->AddHandler("AliReconstructor", detName,
1064 recName, detName, recName + "()");
c757bafd 1065 }
b8cd5251 1066 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1067 }
1068 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1069 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1070 }
b8cd5251 1071 if (reconstructor) {
1072 TObject* obj = fOptions.FindObject(detName.Data());
1073 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1074 reconstructor->Init(fRunLoader);
b8cd5251 1075 fReconstructor[iDet] = reconstructor;
1076 }
1077
f08fc9f5 1078 // get or create the loader
1079 if (detName != "HLT") {
1080 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1081 if (!fLoader[iDet]) {
1082 AliConfig::Instance()
1083 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1084 detName, detName);
1085 // first check if a plugin is defined for the loader
1086 TPluginHandler* pluginHandler =
1087 pluginManager->FindHandler("AliLoader", detName);
1088 // if not, add a plugin for it
1089 if (!pluginHandler) {
1090 TString loaderName = "Ali" + detName + "Loader";
1091 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1092 pluginManager->AddHandler("AliLoader", detName,
1093 loaderName, detName + "base",
1094 loaderName + "(const char*, TFolder*)");
1095 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1096 }
1097 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1098 fLoader[iDet] =
1099 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1100 fRunLoader->GetEventFolder());
1101 }
1102 if (!fLoader[iDet]) { // use default loader
1103 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1104 }
1105 if (!fLoader[iDet]) {
1106 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1107 if (fStopOnError) return NULL;
f08fc9f5 1108 } else {
1109 fRunLoader->AddLoader(fLoader[iDet]);
1110 fRunLoader->CdGAFile();
1111 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1112 fRunLoader->Write(0, TObject::kOverwrite);
1113 }
1114 }
1115 }
1116
b8cd5251 1117 return reconstructor;
c757bafd 1118}
1119
1120//_____________________________________________________________________________
2257f27e 1121Bool_t AliReconstruction::CreateVertexer()
1122{
1123// create the vertexer
1124
b8cd5251 1125 fVertexer = NULL;
1126 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1127 if (itsReconstructor) {
b8cd5251 1128 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1129 }
b8cd5251 1130 if (!fVertexer) {
815c2b38 1131 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1132 if (fStopOnError) return kFALSE;
1133 }
1134
1135 return kTRUE;
1136}
1137
1138//_____________________________________________________________________________
b8cd5251 1139Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1140{
f08fc9f5 1141// create the trackers
24f7a148 1142
b8cd5251 1143 TString detStr = detectors;
1144 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1145 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1146 AliReconstructor* reconstructor = GetReconstructor(iDet);
1147 if (!reconstructor) continue;
1148 TString detName = fgkDetectorName[iDet];
1f46a9ae 1149 if (detName == "HLT") {
1150 fRunHLTTracking = kTRUE;
1151 continue;
1152 }
f08fc9f5 1153
1154 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1155 if (!fTracker[iDet] && (iDet < 7)) {
1156 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1157 if (fStopOnError) return kFALSE;
1158 }
1159 }
1160
24f7a148 1161 return kTRUE;
1162}
1163
1164//_____________________________________________________________________________
b26c3770 1165void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1166{
1167// delete trackers and the run loader and close and delete the file
1168
b8cd5251 1169 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1170 delete fReconstructor[iDet];
1171 fReconstructor[iDet] = NULL;
1172 fLoader[iDet] = NULL;
1173 delete fTracker[iDet];
1174 fTracker[iDet] = NULL;
1175 }
1176 delete fVertexer;
1177 fVertexer = NULL;
e583c30d 1178
1179 delete fRunLoader;
1180 fRunLoader = NULL;
b649205a 1181 delete fRawReader;
1182 fRawReader = NULL;
e583c30d 1183
1184 if (file) {
1185 file->Close();
1186 delete file;
1187 }
b26c3770 1188
1189 if (fileOld) {
1190 fileOld->Close();
1191 delete fileOld;
1192 gSystem->Unlink("AliESDs.old.root");
1193 }
e583c30d 1194}
24f7a148 1195
1196
1197//_____________________________________________________________________________
1198Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1199{
1200// read the ESD event from a file
1201
1202 if (!esd) return kFALSE;
1203 char fileName[256];
1204 sprintf(fileName, "ESD_%d.%d_%s.root",
1205 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1206 if (gSystem->AccessPathName(fileName)) return kFALSE;
1207
f3a97c86 1208 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1209 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1210 TFile* file = TFile::Open(fileName);
1211 if (!file || !file->IsOpen()) {
815c2b38 1212 AliError(Form("opening %s failed", fileName));
24f7a148 1213 delete file;
1214 return kFALSE;
1215 }
1216
1217 gROOT->cd();
1218 delete esd;
1219 esd = (AliESD*) file->Get("ESD");
1220 file->Close();
1221 delete file;
1222 return kTRUE;
1223}
1224
1225//_____________________________________________________________________________
1226void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1227{
1228// write the ESD event to a file
1229
1230 if (!esd) return;
1231 char fileName[256];
1232 sprintf(fileName, "ESD_%d.%d_%s.root",
1233 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1234
815c2b38 1235 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1236 TFile* file = TFile::Open(fileName, "recreate");
1237 if (!file || !file->IsOpen()) {
815c2b38 1238 AliError(Form("opening %s failed", fileName));
24f7a148 1239 } else {
1240 esd->Write("ESD");
1241 file->Close();
1242 }
1243 delete file;
1244}
f3a97c86 1245
1246
1247
1248
1249//_____________________________________________________________________________
1250void AliReconstruction::CreateTag(TFile* file)
1251{
2bdb9d38 1252 /////////////
1253 //muon code//
1254 ////////////
56982dd1 1255 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1256 //Variables
56982dd1 1257 Double_t fX,fY,fZ ;
1258 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1259 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1260 Int_t fCharge;
1261 TLorentzVector fEPvector;
1262
1263 Float_t fZVertexCut = 10.0;
1264 Float_t fRhoVertexCut = 2.0;
1265
1266 Float_t fLowPtCut = 1.0;
1267 Float_t fHighPtCut = 3.0;
1268 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1269 ////////////
1270
1271 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1272
089bf903 1273 // Creates the tags for all the events in a given ESD file
4302e20f 1274 Int_t ntrack;
089bf903 1275 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1276 Int_t nPos, nNeg, nNeutr;
1277 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1278 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1279 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1280 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1281 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1282 Int_t fVertexflag;
1283 TString fVertexName;
4302e20f 1284
f3a97c86 1285 AliRunTag *tag = new AliRunTag();
f3a97c86 1286 AliEventTag *evTag = new AliEventTag();
1287 TTree ttag("T","A Tree with event tags");
2bdb9d38 1288 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1289 btag->SetCompressionLevel(9);
2bdb9d38 1290
f3a97c86 1291 AliInfo(Form("Creating the tags......."));
1292
1293 if (!file || !file->IsOpen()) {
1294 AliError(Form("opening failed"));
1295 delete file;
1296 return ;
2bdb9d38 1297 }
1298 Int_t firstEvent = 0,lastEvent = 0;
f3a97c86 1299 TTree *t = (TTree*) file->Get("esdTree");
1300 TBranch * b = t->GetBranch("ESD");
1301 AliESD *esd = 0;
1302 b->SetAddress(&esd);
2bdb9d38 1303
f3a97c86 1304 tag->SetRunId(esd->GetRunNumber());
2bdb9d38 1305
089bf903 1306 Int_t iNumberOfEvents = b->GetEntries();
2bdb9d38 1307 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1308 ntrack = 0;
1309 nPos = 0;
1310 nNeg = 0;
1311 nNeutr =0;
1312 nK0s = 0;
1313 nNeutrons = 0;
1314 nPi0s = 0;
1315 nGamas = 0;
1316 nProtons = 0;
1317 nKaons = 0;
1318 nPions = 0;
1319 nMuons = 0;
1320 nElectrons = 0;
1321 nCh1GeV = 0;
1322 nCh3GeV = 0;
1323 nCh10GeV = 0;
1324 nMu1GeV = 0;
1325 nMu3GeV = 0;
1326 nMu10GeV = 0;
1327 nEl1GeV = 0;
1328 nEl3GeV = 0;
1329 nEl10GeV = 0;
1330 maxPt = .0;
1331 meanPt = .0;
1332 totalP = .0;
d1a50cb5 1333 fVertexflag = 1;
1334
2bdb9d38 1335 b->GetEntry(iEventNumber);
1336 const AliESDVertex * vertexIn = esd->GetVertex();
d1a50cb5 1337 fVertexName = vertexIn->GetName();
1338 if(fVertexName == "default") fVertexflag = 0;
2bdb9d38 1339 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1340 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1341 UInt_t status = esdTrack->GetStatus();
f3a97c86 1342
2bdb9d38 1343 //select only tracks with ITS refit
1344 if ((status&AliESDtrack::kITSrefit)==0) continue;
1345 //select only tracks with TPC refit
1346 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1347
2bdb9d38 1348 //select only tracks with the "combined PID"
1349 if ((status&AliESDtrack::kESDpid)==0) continue;
1350 Double_t p[3];
1351 esdTrack->GetPxPyPz(p);
1352 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1353 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1354 totalP += momentum;
1355 meanPt += fPt;
1356 if(fPt > maxPt) maxPt = fPt;
4302e20f 1357
2bdb9d38 1358 if(esdTrack->GetSign() > 0) {
1359 nPos++;
56982dd1 1360 if(fPt > fLowPtCut) nCh1GeV++;
1361 if(fPt > fHighPtCut) nCh3GeV++;
1362 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1363 }
1364 if(esdTrack->GetSign() < 0) {
1365 nNeg++;
56982dd1 1366 if(fPt > fLowPtCut) nCh1GeV++;
1367 if(fPt > fHighPtCut) nCh3GeV++;
1368 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1369 }
1370 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1371
2bdb9d38 1372 //PID
1373 Double_t prob[5];
1374 esdTrack->GetESDpid(prob);
4302e20f 1375
2bdb9d38 1376 Double_t rcc = 0.0;
1377 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1378 if(rcc == 0.0) continue;
1379 //Bayes' formula
1380 Double_t w[5];
1381 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1382
2bdb9d38 1383 //protons
1384 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1385 //kaons
1386 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1387 //pions
1388 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1389 //electrons
1390 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1391 nElectrons++;
56982dd1 1392 if(fPt > fLowPtCut) nEl1GeV++;
1393 if(fPt > fHighPtCut) nEl3GeV++;
1394 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1395 }
1396 ntrack++;
1397 }//track loop
1398
1399 /////////////
1400 //muon code//
1401 ////////////
1402 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1403 // loop over all reconstructed tracks (also first track of combination)
1404 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1405 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1406 if (muonTrack == 0x0) continue;
4302e20f 1407
2bdb9d38 1408 // Coordinates at vertex
56982dd1 1409 fZ = muonTrack->GetZ();
1410 fY = muonTrack->GetBendingCoor();
1411 fX = muonTrack->GetNonBendingCoor();
4302e20f 1412
56982dd1 1413 fThetaX = muonTrack->GetThetaX();
1414 fThetaY = muonTrack->GetThetaY();
4302e20f 1415
56982dd1 1416 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1417 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1418 fPxRec = fPzRec * TMath::Tan(fThetaX);
1419 fPyRec = fPzRec * TMath::Tan(fThetaY);
1420 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1421
2bdb9d38 1422 //ChiSquare of the track if needed
56982dd1 1423 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1424 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1425 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1426
2bdb9d38 1427 // total number of muons inside a vertex cut
56982dd1 1428 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1429 nMuons++;
56982dd1 1430 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1431 nMu1GeV++;
56982dd1 1432 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1433 nMu3GeV++;
56982dd1 1434 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1435 nMu10GeV++;
1436 }
1437 }
1438 }
1439 }
1440 }//muon track loop
1441
1442 // Fill the event tags
1443 if(ntrack != 0)
1444 meanPt = meanPt/ntrack;
1445
1446 evTag->SetEventId(iEventNumber+1);
1447 evTag->SetVertexX(vertexIn->GetXv());
1448 evTag->SetVertexY(vertexIn->GetYv());
1449 evTag->SetVertexZ(vertexIn->GetZv());
d1a50cb5 1450 evTag->SetVertexZError(vertexIn->GetZRes());
1451 evTag->SetVertexFlag(fVertexflag);
1452
2bdb9d38 1453 evTag->SetT0VertexZ(esd->GetT0zVertex());
1454
1455 evTag->SetTrigger(esd->GetTrigger());
1456
32a5cab4 1457 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1458 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1459 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1460 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1461 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1462 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1463
1464
1465 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1466 evTag->SetNumOfPosTracks(nPos);
1467 evTag->SetNumOfNegTracks(nNeg);
1468 evTag->SetNumOfNeutrTracks(nNeutr);
1469
1470 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1471 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1472 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1473 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1474
1475 evTag->SetNumOfProtons(nProtons);
1476 evTag->SetNumOfKaons(nKaons);
1477 evTag->SetNumOfPions(nPions);
1478 evTag->SetNumOfMuons(nMuons);
1479 evTag->SetNumOfElectrons(nElectrons);
1480 evTag->SetNumOfPhotons(nGamas);
1481 evTag->SetNumOfPi0s(nPi0s);
1482 evTag->SetNumOfNeutrons(nNeutrons);
1483 evTag->SetNumOfKaon0s(nK0s);
1484
1485 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1486 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1487 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1488 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1489 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1490 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1491 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1492 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1493 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1494
85c60a8e 1495 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1496 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1497
1498 evTag->SetTotalMomentum(totalP);
1499 evTag->SetMeanPt(meanPt);
1500 evTag->SetMaxPt(maxPt);
1501
1502 tag->AddEventTag(*evTag);
1503 }
089bf903 1504 lastEvent = iNumberOfEvents;
f3a97c86 1505
1506 ttag.Fill();
1507 tag->Clear();
1508
1509 char fileName[256];
1510 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1511 tag->GetRunId(),firstEvent,lastEvent );
1512 AliInfo(Form("writing tags to file %s", fileName));
1513 AliDebug(1, Form("writing tags to file %s", fileName));
1514
1515 TFile* ftag = TFile::Open(fileName, "recreate");
1516 ftag->cd();
1517 ttag.Write();
1518 ftag->Close();
1519 file->cd();
1520 delete tag;
f3a97c86 1521 delete evTag;
1522}
1523
98937d93 1524void AliReconstruction::WriteAlignmentData(AliESD* esd)
1525{
1526 // Write space-points which are then used in the alignment procedures
1527 // For the moment only ITS, TRD and TPC
1528
1529 // Load TOF clusters
d528ee75 1530 if (fTracker[3]){
1531 fLoader[3]->LoadRecPoints("read");
1532 TTree* tree = fLoader[3]->TreeR();
1533 if (!tree) {
1534 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1535 return;
1536 }
1537 fTracker[3]->LoadClusters(tree);
98937d93 1538 }
98937d93 1539 Int_t ntracks = esd->GetNumberOfTracks();
1540 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1541 {
1542 AliESDtrack *track = esd->GetTrack(itrack);
1543 Int_t nsp = 0;
ef7253ac 1544 Int_t idx[200];
98937d93 1545 for (Int_t iDet = 3; iDet >= 0; iDet--)
1546 nsp += track->GetNcls(iDet);
1547 if (nsp) {
1548 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1549 track->SetTrackPointArray(sp);
1550 Int_t isptrack = 0;
1551 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1552 AliTracker *tracker = fTracker[iDet];
1553 if (!tracker) continue;
1554 Int_t nspdet = track->GetNcls(iDet);
98937d93 1555 if (nspdet <= 0) continue;
1556 track->GetClusters(iDet,idx);
1557 AliTrackPoint p;
1558 Int_t isp = 0;
1559 Int_t isp2 = 0;
1560 while (isp < nspdet) {
1561 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 1562 const Int_t kNTPCmax = 159;
1563 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 1564 if (!isvalid) continue;
1565 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1566 }
98937d93 1567 }
1568 }
1569 }
d528ee75 1570 if (fTracker[3]){
1571 fTracker[3]->UnloadClusters();
1572 fLoader[3]->UnloadRecPoints();
1573 }
98937d93 1574}