]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliReconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Using AliPID (Yu.Belikov)
[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
c757bafd 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
596a855f 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
b26c3770 540//_____________________________________________________________________________
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
596a855f 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();
2257f27e 639 }
640 else {
641 vertex->SetTruePos(vtxPos); // store also the vertex from MC
642 }
643
644 } else {
815c2b38 645 AliInfo("getting the primary vertex from MC");
2257f27e 646 vertex = new AliESDVertex(vtxPos, vtxErr);
647 }
648
649 if (vertex) {
650 vertex->GetXYZ(vtxPos);
651 vertex->GetSigmaXYZ(vtxErr);
652 } else {
815c2b38 653 AliWarning("no vertex reconstructed");
2257f27e 654 vertex = new AliESDVertex(vtxPos, vtxErr);
655 }
656 esd->SetVertex(vertex);
b8cd5251 657 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
658 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
659 }
2257f27e 660 delete vertex;
661
5f8272e1 662 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
663 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 664
665 return kTRUE;
666}
667
1f46a9ae 668//_____________________________________________________________________________
669Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
670{
671// run the HLT barrel tracking
672
673 TStopwatch stopwatch;
674 stopwatch.Start();
675
676 if (!fRunLoader) {
677 AliError("Missing runLoader!");
678 return kFALSE;
679 }
680
681 AliInfo("running HLT tracking");
682
683 // Get a pointer to the HLT reconstructor
684 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
685 if (!reconstructor) return kFALSE;
686
687 // TPC + ITS
688 for (Int_t iDet = 1; iDet >= 0; iDet--) {
689 TString detName = fgkDetectorName[iDet];
690 AliDebug(1, Form("%s HLT tracking", detName.Data()));
691 reconstructor->SetOption(detName.Data());
692 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
693 if (!tracker) {
694 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
695 if (fStopOnError) return kFALSE;
9dcc06e1 696 continue;
1f46a9ae 697 }
698 Double_t vtxPos[3];
699 Double_t vtxErr[3]={0.005,0.005,0.010};
700 const AliESDVertex *vertex = esd->GetVertex();
701 vertex->GetXYZ(vtxPos);
702 tracker->SetVertex(vtxPos,vtxErr);
703 if(iDet != 1) {
704 fLoader[iDet]->LoadRecPoints("read");
705 TTree* tree = fLoader[iDet]->TreeR();
706 if (!tree) {
707 AliError(Form("Can't get the %s cluster tree", detName.Data()));
708 return kFALSE;
709 }
710 tracker->LoadClusters(tree);
711 }
712 if (tracker->Clusters2Tracks(esd) != 0) {
713 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
714 return kFALSE;
715 }
716 if(iDet != 1) {
717 tracker->UnloadClusters();
718 }
719 delete tracker;
720 }
721
5f8272e1 722 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
723 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 724
725 return kTRUE;
726}
727
2257f27e 728//_____________________________________________________________________________
729Bool_t AliReconstruction::RunTracking(AliESD*& esd)
730{
731// run the barrel tracking
732
733 TStopwatch stopwatch;
734 stopwatch.Start();
24f7a148 735
815c2b38 736 AliInfo("running tracking");
596a855f 737
b8cd5251 738 // pass 1: TPC + ITS inwards
739 for (Int_t iDet = 1; iDet >= 0; iDet--) {
740 if (!fTracker[iDet]) continue;
741 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 742
b8cd5251 743 // load clusters
744 fLoader[iDet]->LoadRecPoints("read");
745 TTree* tree = fLoader[iDet]->TreeR();
746 if (!tree) {
747 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 748 return kFALSE;
749 }
b8cd5251 750 fTracker[iDet]->LoadClusters(tree);
751
752 // run tracking
753 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
754 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 755 return kFALSE;
756 }
b8cd5251 757 if (fCheckPointLevel > 1) {
758 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
759 }
878e1fe1 760 // preliminary PID in TPC needed by the ITS tracker
761 if (iDet == 1) {
762 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 763 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 764 AliESDpid::MakePID(esd);
765 }
b8cd5251 766 }
596a855f 767
b8cd5251 768 // pass 2: ALL backwards
769 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
770 if (!fTracker[iDet]) continue;
771 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
772
773 // load clusters
774 if (iDet > 1) { // all except ITS, TPC
775 TTree* tree = NULL;
7b61cd9c 776 fLoader[iDet]->LoadRecPoints("read");
777 tree = fLoader[iDet]->TreeR();
b8cd5251 778 if (!tree) {
779 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 780 return kFALSE;
781 }
b8cd5251 782 fTracker[iDet]->LoadClusters(tree);
783 }
24f7a148 784
b8cd5251 785 // run tracking
786 if (fTracker[iDet]->PropagateBack(esd) != 0) {
787 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
788 return kFALSE;
789 }
790 if (fCheckPointLevel > 1) {
791 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
792 }
24f7a148 793
b8cd5251 794 // unload clusters
795 if (iDet > 2) { // all except ITS, TPC, TRD
796 fTracker[iDet]->UnloadClusters();
7b61cd9c 797 fLoader[iDet]->UnloadRecPoints();
b8cd5251 798 }
8f37df88 799 // updated PID in TPC needed by the ITS tracker -MI
800 if (iDet == 1) {
801 GetReconstructor(1)->FillESD(fRunLoader, esd);
802 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
803 AliESDpid::MakePID(esd);
804 }
b8cd5251 805 }
596a855f 806
98937d93 807 // write space-points to the ESD in case alignment data output
808 // is switched on
809 if (fWriteAlignmentData)
810 WriteAlignmentData(esd);
811
b8cd5251 812 // pass 3: TRD + TPC + ITS refit inwards
813 for (Int_t iDet = 2; iDet >= 0; iDet--) {
814 if (!fTracker[iDet]) continue;
815 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 816
b8cd5251 817 // run tracking
818 if (fTracker[iDet]->RefitInward(esd) != 0) {
819 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
820 return kFALSE;
821 }
822 if (fCheckPointLevel > 1) {
823 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
824 }
596a855f 825
b8cd5251 826 // unload clusters
827 fTracker[iDet]->UnloadClusters();
828 fLoader[iDet]->UnloadRecPoints();
829 }
ff8bb5ae 830 //
831 // Propagate track to the vertex - if not done by ITS
832 //
833 Int_t ntracks = esd->GetNumberOfTracks();
834 for (Int_t itrack=0; itrack<ntracks; itrack++){
835 const Double_t kRadius = 3; // beam pipe radius
836 const Double_t kMaxStep = 5; // max step
837 const Double_t kMaxD = 123456; // max distance to prim vertex
838 Double_t fieldZ = AliTracker::GetBz(); //
839 AliESDtrack * track = esd->GetTrack(itrack);
840 if (!track) continue;
841 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
842 track->PropagateTo(kRadius, track->GetMass(),kMaxStep,kTRUE);
843 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
844 }
845
5f8272e1 846 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
847 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 848
596a855f 849 return kTRUE;
850}
851
852//_____________________________________________________________________________
24f7a148 853Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 854{
855// fill the event summary data
856
030b532d 857 TStopwatch stopwatch;
858 stopwatch.Start();
815c2b38 859 AliInfo("filling ESD");
030b532d 860
596a855f 861 TString detStr = detectors;
b8cd5251 862 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
863 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
864 AliReconstructor* reconstructor = GetReconstructor(iDet);
865 if (!reconstructor) continue;
866
867 if (!ReadESD(esd, fgkDetectorName[iDet])) {
868 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 869 TTree* clustersTree = NULL;
870 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
871 fLoader[iDet]->LoadRecPoints("read");
872 clustersTree = fLoader[iDet]->TreeR();
873 if (!clustersTree) {
874 AliError(Form("Can't get the %s clusters tree",
875 fgkDetectorName[iDet]));
876 if (fStopOnError) return kFALSE;
877 }
878 }
879 if (fRawReader && !reconstructor->HasDigitConversion()) {
880 reconstructor->FillESD(fRawReader, clustersTree, esd);
881 } else {
882 TTree* digitsTree = NULL;
883 if (fLoader[iDet]) {
884 fLoader[iDet]->LoadDigits("read");
885 digitsTree = fLoader[iDet]->TreeD();
886 if (!digitsTree) {
887 AliError(Form("Can't get the %s digits tree",
888 fgkDetectorName[iDet]));
889 if (fStopOnError) return kFALSE;
890 }
891 }
892 reconstructor->FillESD(digitsTree, clustersTree, esd);
893 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
894 }
895 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
896 fLoader[iDet]->UnloadRecPoints();
897 }
898
b8cd5251 899 if (fRawReader) {
900 reconstructor->FillESD(fRunLoader, fRawReader, esd);
901 } else {
902 reconstructor->FillESD(fRunLoader, esd);
24f7a148 903 }
b8cd5251 904 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 905 }
906 }
907
908 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 909 AliError(Form("the following detectors were not found: %s",
910 detStr.Data()));
596a855f 911 if (fStopOnError) return kFALSE;
912 }
913
5f8272e1 914 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
915 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 916
596a855f 917 return kTRUE;
918}
919
920
921//_____________________________________________________________________________
922Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
923{
924// check whether detName is contained in detectors
925// if yes, it is removed from detectors
926
927 // check if all detectors are selected
928 if ((detectors.CompareTo("ALL") == 0) ||
929 detectors.BeginsWith("ALL ") ||
930 detectors.EndsWith(" ALL") ||
931 detectors.Contains(" ALL ")) {
932 detectors = "ALL";
933 return kTRUE;
934 }
935
936 // search for the given detector
937 Bool_t result = kFALSE;
938 if ((detectors.CompareTo(detName) == 0) ||
939 detectors.BeginsWith(detName+" ") ||
940 detectors.EndsWith(" "+detName) ||
941 detectors.Contains(" "+detName+" ")) {
942 detectors.ReplaceAll(detName, "");
943 result = kTRUE;
944 }
945
946 // clean up the detectors string
947 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
948 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
949 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
950
951 return result;
952}
e583c30d 953
f08fc9f5 954//_____________________________________________________________________________
955Bool_t AliReconstruction::InitRunLoader()
956{
957// get or create the run loader
958
959 if (gAlice) delete gAlice;
960 gAlice = NULL;
961
b26c3770 962 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
963 // load all base libraries to get the loader classes
964 TString libs = gSystem->GetLibraries();
965 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
966 TString detName = fgkDetectorName[iDet];
967 if (detName == "HLT") continue;
968 if (libs.Contains("lib" + detName + "base.so")) continue;
969 gSystem->Load("lib" + detName + "base.so");
970 }
f08fc9f5 971 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
972 if (!fRunLoader) {
973 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
974 CleanUp();
975 return kFALSE;
976 }
b26c3770 977 fRunLoader->CdGAFile();
978 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
979 if (fRunLoader->LoadgAlice() == 0) {
980 gAlice = fRunLoader->GetAliRun();
c84a5e9e 981 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 982 }
f08fc9f5 983 }
984 if (!gAlice && !fRawReader) {
985 AliError(Form("no gAlice object found in file %s",
986 fGAliceFileName.Data()));
987 CleanUp();
988 return kFALSE;
989 }
990
991 } else { // galice.root does not exist
992 if (!fRawReader) {
993 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
994 CleanUp();
995 return kFALSE;
996 }
997 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
998 AliConfig::GetDefaultEventFolderName(),
999 "recreate");
1000 if (!fRunLoader) {
1001 AliError(Form("could not create run loader in file %s",
1002 fGAliceFileName.Data()));
1003 CleanUp();
1004 return kFALSE;
1005 }
1006 fRunLoader->MakeTree("E");
1007 Int_t iEvent = 0;
1008 while (fRawReader->NextEvent()) {
1009 fRunLoader->SetEventNumber(iEvent);
1010 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1011 iEvent, iEvent);
1012 fRunLoader->MakeTree("H");
1013 fRunLoader->TreeE()->Fill();
1014 iEvent++;
1015 }
1016 fRawReader->RewindEvents();
1017 fRunLoader->WriteHeader("OVERWRITE");
1018 fRunLoader->CdGAFile();
1019 fRunLoader->Write(0, TObject::kOverwrite);
1020// AliTracker::SetFieldMap(???);
1021 }
1022
1023 return kTRUE;
1024}
1025
c757bafd 1026//_____________________________________________________________________________
b8cd5251 1027AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1028{
f08fc9f5 1029// get the reconstructor object and the loader for a detector
c757bafd 1030
b8cd5251 1031 if (fReconstructor[iDet]) return fReconstructor[iDet];
1032
1033 // load the reconstructor object
1034 TPluginManager* pluginManager = gROOT->GetPluginManager();
1035 TString detName = fgkDetectorName[iDet];
1036 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1037 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1038
1039 if (detName == "HLT") {
1040 if (!gROOT->GetClass("AliLevel3")) {
1041 gSystem->Load("libAliL3Src.so");
1042 gSystem->Load("libAliL3Misc.so");
1043 gSystem->Load("libAliL3Hough.so");
1044 gSystem->Load("libAliL3Comp.so");
1045 }
1046 }
1047
1048 AliReconstructor* reconstructor = NULL;
1049 // first check if a plugin is defined for the reconstructor
1050 TPluginHandler* pluginHandler =
1051 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1052 // if not, add a plugin for it
1053 if (!pluginHandler) {
b8cd5251 1054 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1055 TString libs = gSystem->GetLibraries();
1056 if (libs.Contains("lib" + detName + "base.so") ||
1057 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1058 pluginManager->AddHandler("AliReconstructor", detName,
1059 recName, detName + "rec", recName + "()");
1060 } else {
1061 pluginManager->AddHandler("AliReconstructor", detName,
1062 recName, detName, recName + "()");
c757bafd 1063 }
b8cd5251 1064 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1065 }
1066 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1067 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1068 }
b8cd5251 1069 if (reconstructor) {
1070 TObject* obj = fOptions.FindObject(detName.Data());
1071 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1072 reconstructor->Init(fRunLoader);
b8cd5251 1073 fReconstructor[iDet] = reconstructor;
1074 }
1075
f08fc9f5 1076 // get or create the loader
1077 if (detName != "HLT") {
1078 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1079 if (!fLoader[iDet]) {
1080 AliConfig::Instance()
1081 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1082 detName, detName);
1083 // first check if a plugin is defined for the loader
1084 TPluginHandler* pluginHandler =
1085 pluginManager->FindHandler("AliLoader", detName);
1086 // if not, add a plugin for it
1087 if (!pluginHandler) {
1088 TString loaderName = "Ali" + detName + "Loader";
1089 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1090 pluginManager->AddHandler("AliLoader", detName,
1091 loaderName, detName + "base",
1092 loaderName + "(const char*, TFolder*)");
1093 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1094 }
1095 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1096 fLoader[iDet] =
1097 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1098 fRunLoader->GetEventFolder());
1099 }
1100 if (!fLoader[iDet]) { // use default loader
1101 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1102 }
1103 if (!fLoader[iDet]) {
1104 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1105 if (fStopOnError) return NULL;
f08fc9f5 1106 } else {
1107 fRunLoader->AddLoader(fLoader[iDet]);
1108 fRunLoader->CdGAFile();
1109 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1110 fRunLoader->Write(0, TObject::kOverwrite);
1111 }
1112 }
1113 }
1114
b8cd5251 1115 return reconstructor;
c757bafd 1116}
1117
2257f27e 1118//_____________________________________________________________________________
1119Bool_t AliReconstruction::CreateVertexer()
1120{
1121// create the vertexer
1122
b8cd5251 1123 fVertexer = NULL;
1124 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1125 if (itsReconstructor) {
b8cd5251 1126 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1127 }
b8cd5251 1128 if (!fVertexer) {
815c2b38 1129 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1130 if (fStopOnError) return kFALSE;
1131 }
1132
1133 return kTRUE;
1134}
1135
24f7a148 1136//_____________________________________________________________________________
b8cd5251 1137Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1138{
f08fc9f5 1139// create the trackers
24f7a148 1140
b8cd5251 1141 TString detStr = detectors;
1142 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1143 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1144 AliReconstructor* reconstructor = GetReconstructor(iDet);
1145 if (!reconstructor) continue;
1146 TString detName = fgkDetectorName[iDet];
1f46a9ae 1147 if (detName == "HLT") {
1148 fRunHLTTracking = kTRUE;
1149 continue;
1150 }
f08fc9f5 1151
1152 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1153 if (!fTracker[iDet] && (iDet < 7)) {
1154 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1155 if (fStopOnError) return kFALSE;
1156 }
1157 }
1158
24f7a148 1159 return kTRUE;
1160}
1161
e583c30d 1162//_____________________________________________________________________________
b26c3770 1163void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1164{
1165// delete trackers and the run loader and close and delete the file
1166
b8cd5251 1167 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1168 delete fReconstructor[iDet];
1169 fReconstructor[iDet] = NULL;
1170 fLoader[iDet] = NULL;
1171 delete fTracker[iDet];
1172 fTracker[iDet] = NULL;
1173 }
1174 delete fVertexer;
1175 fVertexer = NULL;
e583c30d 1176
1177 delete fRunLoader;
1178 fRunLoader = NULL;
b649205a 1179 delete fRawReader;
1180 fRawReader = NULL;
e583c30d 1181
1182 if (file) {
1183 file->Close();
1184 delete file;
1185 }
b26c3770 1186
1187 if (fileOld) {
1188 fileOld->Close();
1189 delete fileOld;
1190 gSystem->Unlink("AliESDs.old.root");
1191 }
e583c30d 1192}
24f7a148 1193
1194
1195//_____________________________________________________________________________
1196Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1197{
1198// read the ESD event from a file
1199
1200 if (!esd) return kFALSE;
1201 char fileName[256];
1202 sprintf(fileName, "ESD_%d.%d_%s.root",
1203 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1204 if (gSystem->AccessPathName(fileName)) return kFALSE;
1205
f3a97c86 1206 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1207 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1208 TFile* file = TFile::Open(fileName);
1209 if (!file || !file->IsOpen()) {
815c2b38 1210 AliError(Form("opening %s failed", fileName));
24f7a148 1211 delete file;
1212 return kFALSE;
1213 }
1214
1215 gROOT->cd();
1216 delete esd;
1217 esd = (AliESD*) file->Get("ESD");
1218 file->Close();
1219 delete file;
1220 return kTRUE;
1221}
1222
1223//_____________________________________________________________________________
1224void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1225{
1226// write the ESD event to a file
1227
1228 if (!esd) return;
1229 char fileName[256];
1230 sprintf(fileName, "ESD_%d.%d_%s.root",
1231 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1232
815c2b38 1233 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1234 TFile* file = TFile::Open(fileName, "recreate");
1235 if (!file || !file->IsOpen()) {
815c2b38 1236 AliError(Form("opening %s failed", fileName));
24f7a148 1237 } else {
1238 esd->Write("ESD");
1239 file->Close();
1240 }
1241 delete file;
1242}
f3a97c86 1243
1244
1245
1246
1247//_____________________________________________________________________________
1248void AliReconstruction::CreateTag(TFile* file)
1249{
2bdb9d38 1250 /////////////
1251 //muon code//
1252 ////////////
56982dd1 1253 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1254 //Variables
56982dd1 1255 Double_t fX,fY,fZ ;
1256 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1257 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1258 Int_t fCharge;
1259 TLorentzVector fEPvector;
1260
1261 Float_t fZVertexCut = 10.0;
1262 Float_t fRhoVertexCut = 2.0;
1263
1264 Float_t fLowPtCut = 1.0;
1265 Float_t fHighPtCut = 3.0;
1266 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1267 ////////////
1268
1269 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1270
089bf903 1271 // Creates the tags for all the events in a given ESD file
4302e20f 1272 Int_t ntrack;
089bf903 1273 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1274 Int_t nPos, nNeg, nNeutr;
1275 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1276 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1277 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1278 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1279 Float_t maxPt = .0, meanPt = .0, totalP = .0;
4302e20f 1280
f3a97c86 1281 AliRunTag *tag = new AliRunTag();
f3a97c86 1282 AliEventTag *evTag = new AliEventTag();
1283 TTree ttag("T","A Tree with event tags");
2bdb9d38 1284 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1285 btag->SetCompressionLevel(9);
2bdb9d38 1286
f3a97c86 1287 AliInfo(Form("Creating the tags......."));
1288
1289 if (!file || !file->IsOpen()) {
1290 AliError(Form("opening failed"));
1291 delete file;
1292 return ;
2bdb9d38 1293 }
1294 Int_t firstEvent = 0,lastEvent = 0;
f3a97c86 1295 TTree *t = (TTree*) file->Get("esdTree");
1296 TBranch * b = t->GetBranch("ESD");
1297 AliESD *esd = 0;
1298 b->SetAddress(&esd);
2bdb9d38 1299
f3a97c86 1300 tag->SetRunId(esd->GetRunNumber());
2bdb9d38 1301
089bf903 1302 Int_t iNumberOfEvents = b->GetEntries();
2bdb9d38 1303 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1304 ntrack = 0;
1305 nPos = 0;
1306 nNeg = 0;
1307 nNeutr =0;
1308 nK0s = 0;
1309 nNeutrons = 0;
1310 nPi0s = 0;
1311 nGamas = 0;
1312 nProtons = 0;
1313 nKaons = 0;
1314 nPions = 0;
1315 nMuons = 0;
1316 nElectrons = 0;
1317 nCh1GeV = 0;
1318 nCh3GeV = 0;
1319 nCh10GeV = 0;
1320 nMu1GeV = 0;
1321 nMu3GeV = 0;
1322 nMu10GeV = 0;
1323 nEl1GeV = 0;
1324 nEl3GeV = 0;
1325 nEl10GeV = 0;
1326 maxPt = .0;
1327 meanPt = .0;
1328 totalP = .0;
1329
1330 b->GetEntry(iEventNumber);
1331 const AliESDVertex * vertexIn = esd->GetVertex();
1332
1333 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1334 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1335 UInt_t status = esdTrack->GetStatus();
f3a97c86 1336
2bdb9d38 1337 //select only tracks with ITS refit
1338 if ((status&AliESDtrack::kITSrefit)==0) continue;
1339 //select only tracks with TPC refit
1340 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1341
2bdb9d38 1342 //select only tracks with the "combined PID"
1343 if ((status&AliESDtrack::kESDpid)==0) continue;
1344 Double_t p[3];
1345 esdTrack->GetPxPyPz(p);
1346 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1347 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1348 totalP += momentum;
1349 meanPt += fPt;
1350 if(fPt > maxPt) maxPt = fPt;
4302e20f 1351
2bdb9d38 1352 if(esdTrack->GetSign() > 0) {
1353 nPos++;
56982dd1 1354 if(fPt > fLowPtCut) nCh1GeV++;
1355 if(fPt > fHighPtCut) nCh3GeV++;
1356 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1357 }
1358 if(esdTrack->GetSign() < 0) {
1359 nNeg++;
56982dd1 1360 if(fPt > fLowPtCut) nCh1GeV++;
1361 if(fPt > fHighPtCut) nCh3GeV++;
1362 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1363 }
1364 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1365
2bdb9d38 1366 //PID
1367 Double_t prob[5];
1368 esdTrack->GetESDpid(prob);
4302e20f 1369
2bdb9d38 1370 Double_t rcc = 0.0;
1371 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1372 if(rcc == 0.0) continue;
1373 //Bayes' formula
1374 Double_t w[5];
1375 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1376
2bdb9d38 1377 //protons
1378 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1379 //kaons
1380 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1381 //pions
1382 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1383 //electrons
1384 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1385 nElectrons++;
56982dd1 1386 if(fPt > fLowPtCut) nEl1GeV++;
1387 if(fPt > fHighPtCut) nEl3GeV++;
1388 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1389 }
1390 ntrack++;
1391 }//track loop
1392
1393 /////////////
1394 //muon code//
1395 ////////////
1396 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1397 // loop over all reconstructed tracks (also first track of combination)
1398 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1399 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1400 if (muonTrack == 0x0) continue;
4302e20f 1401
2bdb9d38 1402 // Coordinates at vertex
56982dd1 1403 fZ = muonTrack->GetZ();
1404 fY = muonTrack->GetBendingCoor();
1405 fX = muonTrack->GetNonBendingCoor();
4302e20f 1406
56982dd1 1407 fThetaX = muonTrack->GetThetaX();
1408 fThetaY = muonTrack->GetThetaY();
4302e20f 1409
56982dd1 1410 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1411 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1412 fPxRec = fPzRec * TMath::Tan(fThetaX);
1413 fPyRec = fPzRec * TMath::Tan(fThetaY);
1414 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1415
2bdb9d38 1416 //ChiSquare of the track if needed
56982dd1 1417 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1418 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1419 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1420
2bdb9d38 1421 // total number of muons inside a vertex cut
56982dd1 1422 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1423 nMuons++;
56982dd1 1424 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1425 nMu1GeV++;
56982dd1 1426 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1427 nMu3GeV++;
56982dd1 1428 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1429 nMu10GeV++;
1430 }
1431 }
1432 }
1433 }
1434 }//muon track loop
1435
1436 // Fill the event tags
1437 if(ntrack != 0)
1438 meanPt = meanPt/ntrack;
1439
1440 evTag->SetEventId(iEventNumber+1);
1441 evTag->SetVertexX(vertexIn->GetXv());
1442 evTag->SetVertexY(vertexIn->GetYv());
1443 evTag->SetVertexZ(vertexIn->GetZv());
1444
1445 evTag->SetT0VertexZ(esd->GetT0zVertex());
1446
1447 evTag->SetTrigger(esd->GetTrigger());
1448
32a5cab4 1449 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1450 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1451 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1452 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1453 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1454 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1455
1456
1457 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1458 evTag->SetNumOfPosTracks(nPos);
1459 evTag->SetNumOfNegTracks(nNeg);
1460 evTag->SetNumOfNeutrTracks(nNeutr);
1461
1462 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1463 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1464 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1465 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1466
1467 evTag->SetNumOfProtons(nProtons);
1468 evTag->SetNumOfKaons(nKaons);
1469 evTag->SetNumOfPions(nPions);
1470 evTag->SetNumOfMuons(nMuons);
1471 evTag->SetNumOfElectrons(nElectrons);
1472 evTag->SetNumOfPhotons(nGamas);
1473 evTag->SetNumOfPi0s(nPi0s);
1474 evTag->SetNumOfNeutrons(nNeutrons);
1475 evTag->SetNumOfKaon0s(nK0s);
1476
1477 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1478 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1479 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1480 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1481 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1482 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1483 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1484 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1485 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1486
85c60a8e 1487 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1488 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1489
1490 evTag->SetTotalMomentum(totalP);
1491 evTag->SetMeanPt(meanPt);
1492 evTag->SetMaxPt(maxPt);
1493
1494 tag->AddEventTag(*evTag);
1495 }
089bf903 1496 lastEvent = iNumberOfEvents;
f3a97c86 1497
1498 ttag.Fill();
1499 tag->Clear();
1500
1501 char fileName[256];
1502 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1503 tag->GetRunId(),firstEvent,lastEvent );
1504 AliInfo(Form("writing tags to file %s", fileName));
1505 AliDebug(1, Form("writing tags to file %s", fileName));
1506
1507 TFile* ftag = TFile::Open(fileName, "recreate");
1508 ftag->cd();
1509 ttag.Write();
1510 ftag->Close();
1511 file->cd();
1512 delete tag;
f3a97c86 1513 delete evTag;
1514}
1515
98937d93 1516void AliReconstruction::WriteAlignmentData(AliESD* esd)
1517{
1518 // Write space-points which are then used in the alignment procedures
1519 // For the moment only ITS, TRD and TPC
1520
1521 // Load TOF clusters
d528ee75 1522 if (fTracker[3]){
1523 fLoader[3]->LoadRecPoints("read");
1524 TTree* tree = fLoader[3]->TreeR();
1525 if (!tree) {
1526 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1527 return;
1528 }
1529 fTracker[3]->LoadClusters(tree);
98937d93 1530 }
98937d93 1531 Int_t ntracks = esd->GetNumberOfTracks();
1532 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1533 {
1534 AliESDtrack *track = esd->GetTrack(itrack);
1535 Int_t nsp = 0;
1536 UInt_t idx[200];
1537 for (Int_t iDet = 3; iDet >= 0; iDet--)
1538 nsp += track->GetNcls(iDet);
1539 if (nsp) {
1540 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1541 track->SetTrackPointArray(sp);
1542 Int_t isptrack = 0;
1543 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1544 AliTracker *tracker = fTracker[iDet];
1545 if (!tracker) continue;
1546 Int_t nspdet = track->GetNcls(iDet);
98937d93 1547 if (nspdet <= 0) continue;
1548 track->GetClusters(iDet,idx);
1549 AliTrackPoint p;
1550 Int_t isp = 0;
1551 Int_t isp2 = 0;
1552 while (isp < nspdet) {
1553 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1554 if (!isvalid) continue;
1555 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1556 }
98937d93 1557 }
1558 }
1559 }
d528ee75 1560 if (fTracker[3]){
1561 fTracker[3]->UnloadClusters();
1562 fLoader[3]->UnloadRecPoints();
1563 }
98937d93 1564}