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