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