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