]> 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
Laurent's mapping modified for trigger chambers
[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"
f3a97c86 132
133#include "AliRunTag.h"
089bf903 134//#include "AliLHCTag.h"
f3a97c86 135#include "AliDetectorTag.h"
136#include "AliEventTag.h"
137
98937d93 138#include "AliTrackPointArray.h"
f3a97c86 139
596a855f 140ClassImp(AliReconstruction)
141
142
c757bafd 143//_____________________________________________________________________________
482070f2 144const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
c757bafd 145
596a855f 146//_____________________________________________________________________________
e583c30d 147AliReconstruction::AliReconstruction(const char* gAliceFilename,
148 const char* name, const char* title) :
149 TNamed(name, title),
150
59697224 151 fRunLocalReconstruction("ALL"),
c84a5e9e 152 fUniformField(kTRUE),
2257f27e 153 fRunVertexFinder(kTRUE),
1f46a9ae 154 fRunHLTTracking(kFALSE),
b8cd5251 155 fRunTracking("ALL"),
e583c30d 156 fFillESD("ALL"),
157 fGAliceFileName(gAliceFilename),
b649205a 158 fInput(""),
b26c3770 159 fFirstEvent(0),
160 fLastEvent(-1),
e583c30d 161 fStopOnError(kFALSE),
24f7a148 162 fCheckPointLevel(0),
b8cd5251 163 fOptions(),
e583c30d 164
165 fRunLoader(NULL),
b649205a 166 fRawReader(NULL),
b8cd5251 167
98937d93 168 fVertexer(NULL),
169
170 fWriteAlignmentData(kFALSE)
596a855f 171{
172// create reconstruction object with default parameters
b8cd5251 173
174 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
175 fReconstructor[iDet] = NULL;
176 fLoader[iDet] = NULL;
177 fTracker[iDet] = NULL;
178 }
e47c4c2e 179 AliPID pid;
3103d196 180 // Import TGeo geometry
181 TGeoManager::Import("geometry.root");
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);
b26c3770 917 }
f08fc9f5 918 }
919 if (!gAlice && !fRawReader) {
920 AliError(Form("no gAlice object found in file %s",
921 fGAliceFileName.Data()));
922 CleanUp();
923 return kFALSE;
924 }
925
926 } else { // galice.root does not exist
927 if (!fRawReader) {
928 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
929 CleanUp();
930 return kFALSE;
931 }
932 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
933 AliConfig::GetDefaultEventFolderName(),
934 "recreate");
935 if (!fRunLoader) {
936 AliError(Form("could not create run loader in file %s",
937 fGAliceFileName.Data()));
938 CleanUp();
939 return kFALSE;
940 }
941 fRunLoader->MakeTree("E");
942 Int_t iEvent = 0;
943 while (fRawReader->NextEvent()) {
944 fRunLoader->SetEventNumber(iEvent);
945 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
946 iEvent, iEvent);
947 fRunLoader->MakeTree("H");
948 fRunLoader->TreeE()->Fill();
949 iEvent++;
950 }
951 fRawReader->RewindEvents();
952 fRunLoader->WriteHeader("OVERWRITE");
953 fRunLoader->CdGAFile();
954 fRunLoader->Write(0, TObject::kOverwrite);
955// AliTracker::SetFieldMap(???);
956 }
957
958 return kTRUE;
959}
960
c757bafd 961//_____________________________________________________________________________
b8cd5251 962AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 963{
f08fc9f5 964// get the reconstructor object and the loader for a detector
c757bafd 965
b8cd5251 966 if (fReconstructor[iDet]) return fReconstructor[iDet];
967
968 // load the reconstructor object
969 TPluginManager* pluginManager = gROOT->GetPluginManager();
970 TString detName = fgkDetectorName[iDet];
971 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 972 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 973
974 if (detName == "HLT") {
975 if (!gROOT->GetClass("AliLevel3")) {
976 gSystem->Load("libAliL3Src.so");
977 gSystem->Load("libAliL3Misc.so");
978 gSystem->Load("libAliL3Hough.so");
979 gSystem->Load("libAliL3Comp.so");
980 }
981 }
982
983 AliReconstructor* reconstructor = NULL;
984 // first check if a plugin is defined for the reconstructor
985 TPluginHandler* pluginHandler =
986 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 987 // if not, add a plugin for it
988 if (!pluginHandler) {
b8cd5251 989 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 990 TString libs = gSystem->GetLibraries();
991 if (libs.Contains("lib" + detName + "base.so") ||
992 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 993 pluginManager->AddHandler("AliReconstructor", detName,
994 recName, detName + "rec", recName + "()");
995 } else {
996 pluginManager->AddHandler("AliReconstructor", detName,
997 recName, detName, recName + "()");
c757bafd 998 }
b8cd5251 999 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1000 }
1001 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1002 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1003 }
b8cd5251 1004 if (reconstructor) {
1005 TObject* obj = fOptions.FindObject(detName.Data());
1006 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1007 reconstructor->Init(fRunLoader);
b8cd5251 1008 fReconstructor[iDet] = reconstructor;
1009 }
1010
f08fc9f5 1011 // get or create the loader
1012 if (detName != "HLT") {
1013 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1014 if (!fLoader[iDet]) {
1015 AliConfig::Instance()
1016 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1017 detName, detName);
1018 // first check if a plugin is defined for the loader
1019 TPluginHandler* pluginHandler =
1020 pluginManager->FindHandler("AliLoader", detName);
1021 // if not, add a plugin for it
1022 if (!pluginHandler) {
1023 TString loaderName = "Ali" + detName + "Loader";
1024 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1025 pluginManager->AddHandler("AliLoader", detName,
1026 loaderName, detName + "base",
1027 loaderName + "(const char*, TFolder*)");
1028 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1029 }
1030 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1031 fLoader[iDet] =
1032 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1033 fRunLoader->GetEventFolder());
1034 }
1035 if (!fLoader[iDet]) { // use default loader
1036 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1037 }
1038 if (!fLoader[iDet]) {
1039 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1040 if (fStopOnError) return NULL;
f08fc9f5 1041 } else {
1042 fRunLoader->AddLoader(fLoader[iDet]);
1043 fRunLoader->CdGAFile();
1044 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1045 fRunLoader->Write(0, TObject::kOverwrite);
1046 }
1047 }
1048 }
1049
b8cd5251 1050 return reconstructor;
c757bafd 1051}
1052
2257f27e 1053//_____________________________________________________________________________
1054Bool_t AliReconstruction::CreateVertexer()
1055{
1056// create the vertexer
1057
b8cd5251 1058 fVertexer = NULL;
1059 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1060 if (itsReconstructor) {
b8cd5251 1061 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1062 }
b8cd5251 1063 if (!fVertexer) {
815c2b38 1064 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1065 if (fStopOnError) return kFALSE;
1066 }
1067
1068 return kTRUE;
1069}
1070
24f7a148 1071//_____________________________________________________________________________
b8cd5251 1072Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1073{
f08fc9f5 1074// create the trackers
24f7a148 1075
b8cd5251 1076 TString detStr = detectors;
1077 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1078 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1079 AliReconstructor* reconstructor = GetReconstructor(iDet);
1080 if (!reconstructor) continue;
1081 TString detName = fgkDetectorName[iDet];
1f46a9ae 1082 if (detName == "HLT") {
1083 fRunHLTTracking = kTRUE;
1084 continue;
1085 }
f08fc9f5 1086
1087 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1088 if (!fTracker[iDet] && (iDet < 7)) {
1089 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1090 if (fStopOnError) return kFALSE;
1091 }
1092 }
1093
24f7a148 1094 return kTRUE;
1095}
1096
e583c30d 1097//_____________________________________________________________________________
b26c3770 1098void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1099{
1100// delete trackers and the run loader and close and delete the file
1101
b8cd5251 1102 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1103 delete fReconstructor[iDet];
1104 fReconstructor[iDet] = NULL;
1105 fLoader[iDet] = NULL;
1106 delete fTracker[iDet];
1107 fTracker[iDet] = NULL;
1108 }
1109 delete fVertexer;
1110 fVertexer = NULL;
e583c30d 1111
1112 delete fRunLoader;
1113 fRunLoader = NULL;
b649205a 1114 delete fRawReader;
1115 fRawReader = NULL;
e583c30d 1116
1117 if (file) {
1118 file->Close();
1119 delete file;
1120 }
b26c3770 1121
1122 if (fileOld) {
1123 fileOld->Close();
1124 delete fileOld;
1125 gSystem->Unlink("AliESDs.old.root");
1126 }
e583c30d 1127}
24f7a148 1128
1129
1130//_____________________________________________________________________________
1131Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1132{
1133// read the ESD event from a file
1134
1135 if (!esd) return kFALSE;
1136 char fileName[256];
1137 sprintf(fileName, "ESD_%d.%d_%s.root",
1138 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1139 if (gSystem->AccessPathName(fileName)) return kFALSE;
1140
f3a97c86 1141 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1142 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1143 TFile* file = TFile::Open(fileName);
1144 if (!file || !file->IsOpen()) {
815c2b38 1145 AliError(Form("opening %s failed", fileName));
24f7a148 1146 delete file;
1147 return kFALSE;
1148 }
1149
1150 gROOT->cd();
1151 delete esd;
1152 esd = (AliESD*) file->Get("ESD");
1153 file->Close();
1154 delete file;
1155 return kTRUE;
1156}
1157
1158//_____________________________________________________________________________
1159void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1160{
1161// write the ESD event to a file
1162
1163 if (!esd) return;
1164 char fileName[256];
1165 sprintf(fileName, "ESD_%d.%d_%s.root",
1166 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1167
815c2b38 1168 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1169 TFile* file = TFile::Open(fileName, "recreate");
1170 if (!file || !file->IsOpen()) {
815c2b38 1171 AliError(Form("opening %s failed", fileName));
24f7a148 1172 } else {
1173 esd->Write("ESD");
1174 file->Close();
1175 }
1176 delete file;
1177}
f3a97c86 1178
1179
1180
1181
1182//_____________________________________________________________________________
1183void AliReconstruction::CreateTag(TFile* file)
1184{
089bf903 1185 // Creates the tags for all the events in a given ESD file
4302e20f 1186 Int_t ntrack;
089bf903 1187 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1188 Int_t nPos, nNeg, nNeutr;
1189 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1190 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1191 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1192 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1193 Float_t maxPt = .0, meanPt = .0, totalP = .0;
4302e20f 1194
f3a97c86 1195 AliRunTag *tag = new AliRunTag();
1196 AliDetectorTag *detTag = new AliDetectorTag();
1197 AliEventTag *evTag = new AliEventTag();
1198 TTree ttag("T","A Tree with event tags");
1199 TBranch * btag = ttag.Branch("AliTAG", "AliRunTag", &tag);
1200 btag->SetCompressionLevel(9);
1201
1202 AliInfo(Form("Creating the tags......."));
1203
1204 if (!file || !file->IsOpen()) {
1205 AliError(Form("opening failed"));
1206 delete file;
1207 return ;
1208 }
1209
1210 TTree *t = (TTree*) file->Get("esdTree");
1211 TBranch * b = t->GetBranch("ESD");
1212 AliESD *esd = 0;
1213 b->SetAddress(&esd);
1214
1215 tag->SetRunId(esd->GetRunNumber());
1216
1217 Int_t firstEvent = 0,lastEvent = 0;
089bf903 1218 Int_t iNumberOfEvents = b->GetEntries();
1219 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++)
f3a97c86 1220 {
4302e20f 1221 ntrack = 0;
089bf903 1222 nPos = 0;
1223 nNeg = 0;
1224 nNeutr =0;
1225 nK0s = 0;
1226 nNeutrons = 0;
1227 nPi0s = 0;
1228 nGamas = 0;
1229 nProtons = 0;
1230 nKaons = 0;
1231 nPions = 0;
1232 nMuons = 0;
1233 nElectrons = 0;
1234 nCh1GeV = 0;
1235 nCh3GeV = 0;
1236 nCh10GeV = 0;
1237 nMu1GeV = 0;
1238 nMu3GeV = 0;
1239 nMu10GeV = 0;
1240 nEl1GeV = 0;
1241 nEl3GeV = 0;
1242 nEl10GeV = 0;
1243 maxPt = .0;
1244 meanPt = .0;
1245 totalP = .0;
1246
1247 b->GetEntry(iEventNumber);
1248 const AliESDVertex * vertexIn = esd->GetVertex();
1249
1250 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++)
4302e20f 1251 {
089bf903 1252 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1253 UInt_t status = esdTrack->GetStatus();
4302e20f 1254
1255 //select only tracks with ITS refit
1256 if ((status&AliESDtrack::kITSrefit)==0) continue;
1257
1258 //select only tracks with TPC refit-->remove extremely high Pt tracks
1259 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1260
1261 //select only tracks with the "combined PID"
1262 if ((status&AliESDtrack::kESDpid)==0) continue;
1263 Double_t p[3];
089bf903 1264 esdTrack->GetPxPyPz(p);
1265 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
4302e20f 1266 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
089bf903 1267 totalP += momentum;
1268 meanPt += fPt;
1269 if(fPt > maxPt)
1270 maxPt = fPt;
4302e20f 1271
089bf903 1272 if(esdTrack->GetSign() > 0)
4302e20f 1273 {
089bf903 1274 nPos++;
4302e20f 1275 if(fPt > 1.0)
089bf903 1276 nCh1GeV++;
4302e20f 1277 if(fPt > 3.0)
089bf903 1278 nCh3GeV++;
4302e20f 1279 if(fPt > 10.0)
089bf903 1280 nCh10GeV++;
4302e20f 1281 }
089bf903 1282 if(esdTrack->GetSign() < 0)
4302e20f 1283 {
089bf903 1284 nNeg++;
4302e20f 1285 if(fPt > 1.0)
089bf903 1286 nCh1GeV++;
4302e20f 1287 if(fPt > 3.0)
089bf903 1288 nCh3GeV++;
4302e20f 1289 if(fPt > 10.0)
089bf903 1290 nCh10GeV++;
4302e20f 1291 }
089bf903 1292 if(esdTrack->GetSign() == 0)
1293 nNeutr++;
4302e20f 1294
1295 //PID
1296 Double_t prob[10];
089bf903 1297 esdTrack->GetESDpid(prob);
4302e20f 1298
1299 //K0s
1300 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 1301 nK0s++;
4302e20f 1302 //neutrons
1303 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 1304 nNeutrons++;
4302e20f 1305 //pi0s
1306 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 1307 nPi0s++;
4302e20f 1308 //gamas
1309 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 1310 nGamas++;
4302e20f 1311 //protons
1312 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 1313 nProtons++;
4302e20f 1314 //kaons
1315 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 1316 nKaons++;
4302e20f 1317 //kaons
1318 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 1319 nPions++;
4302e20f 1320 //muons
1321 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]))
1322 {
089bf903 1323 nMuons++;
4302e20f 1324 if(fPt > 1.0)
089bf903 1325 nMu1GeV++;
4302e20f 1326 if(fPt > 3.0)
089bf903 1327 nMu3GeV++;
4302e20f 1328 if(fPt > 10.0)
089bf903 1329 nMu10GeV++;
4302e20f 1330 }
1331 //electrons
1332 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]))
1333 {
089bf903 1334 nElectrons++;
4302e20f 1335 if(fPt > 1.0)
089bf903 1336 nEl1GeV++;
4302e20f 1337 if(fPt > 3.0)
089bf903 1338 nEl3GeV++;
4302e20f 1339 if(fPt > 10.0)
089bf903 1340 nEl10GeV++;
4302e20f 1341 }
1342
1343
1344
1345 ntrack++;
1346 }//track loop
1347 // Fill the event tags
b085f7b9 1348 if(ntrack != 0)
089bf903 1349 meanPt = meanPt/ntrack;
f3a97c86 1350
089bf903 1351 evTag->SetEventId(iEventNumber+1);
1352 evTag->SetVertexX(vertexIn->GetXv());
1353 evTag->SetVertexY(vertexIn->GetYv());
1354 evTag->SetVertexZ(vertexIn->GetZv());
4302e20f 1355
1356 evTag->SetT0VertexZ(esd->GetT0zVertex());
1357
1358 evTag->SetTrigger(esd->GetTrigger());
1359
1360 evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
1361 evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
1362 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1363 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1364
1365
f3a97c86 1366 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
089bf903 1367 evTag->SetNumOfPosTracks(nPos);
1368 evTag->SetNumOfNegTracks(nNeg);
1369 evTag->SetNumOfNeutrTracks(nNeutr);
4302e20f 1370
f3a97c86 1371 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1372 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
4302e20f 1373 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1374 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1375
089bf903 1376 evTag->SetNumOfProtons(nProtons);
1377 evTag->SetNumOfKaons(nKaons);
1378 evTag->SetNumOfPions(nPions);
1379 evTag->SetNumOfMuons(nMuons);
1380 evTag->SetNumOfElectrons(nElectrons);
1381 evTag->SetNumOfPhotons(nGamas);
1382 evTag->SetNumOfPi0s(nPi0s);
1383 evTag->SetNumOfNeutrons(nNeutrons);
1384 evTag->SetNumOfKaon0s(nK0s);
4302e20f 1385
089bf903 1386 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1387 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1388 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1389 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1390 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1391 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1392 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1393 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1394 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
f3a97c86 1395
1396 evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
1397 evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
4302e20f 1398
089bf903 1399 evTag->SetTotalMomentum(totalP);
1400 evTag->SetMeanPt(meanPt);
1401 evTag->SetMaxPt(maxPt);
4302e20f 1402
cb1645b7 1403 tag->AddEventTag(*evTag);
f3a97c86 1404 }
089bf903 1405 lastEvent = iNumberOfEvents;
f3a97c86 1406
1407 ttag.Fill();
1408 tag->Clear();
1409
1410 char fileName[256];
1411 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1412 tag->GetRunId(),firstEvent,lastEvent );
1413 AliInfo(Form("writing tags to file %s", fileName));
1414 AliDebug(1, Form("writing tags to file %s", fileName));
1415
1416 TFile* ftag = TFile::Open(fileName, "recreate");
1417 ftag->cd();
1418 ttag.Write();
1419 ftag->Close();
1420 file->cd();
1421 delete tag;
1422 delete detTag;
1423 delete evTag;
1424}
1425
98937d93 1426void AliReconstruction::WriteAlignmentData(AliESD* esd)
1427{
1428 // Write space-points which are then used in the alignment procedures
1429 // For the moment only ITS, TRD and TPC
1430
1431 // Load TOF clusters
1432 fLoader[3]->LoadRecPoints("read");
1433 TTree* tree = fLoader[3]->TreeR();
1434 if (!tree) {
1435 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1436 return;
1437 }
1438 fTracker[3]->LoadClusters(tree);
1439
1440 Int_t ntracks = esd->GetNumberOfTracks();
1441 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1442 {
1443 AliESDtrack *track = esd->GetTrack(itrack);
1444 Int_t nsp = 0;
1445 UInt_t idx[200];
1446 for (Int_t iDet = 3; iDet >= 0; iDet--)
1447 nsp += track->GetNcls(iDet);
1448 if (nsp) {
1449 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1450 track->SetTrackPointArray(sp);
1451 Int_t isptrack = 0;
1452 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1453 AliTracker *tracker = fTracker[iDet];
1454 if (!tracker) continue;
1455 Int_t nspdet = track->GetNcls(iDet);
1456 cout<<iDet<<" "<<nspdet<<endl;
1457 if (nspdet <= 0) continue;
1458 track->GetClusters(iDet,idx);
1459 AliTrackPoint p;
1460 Int_t isp = 0;
1461 Int_t isp2 = 0;
1462 while (isp < nspdet) {
1463 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1464 if (!isvalid) continue;
1465 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1466 }
1467 // for (Int_t isp = 0; isp < nspdet; isp++) {
1468 // AliCluster *cl = tracker->GetCluster(idx[isp]);
1469 // UShort_t volid = tracker->GetVolumeID(idx[isp]);
1470 // tracker->GetTrackPoint(idx[isp],p);
1471 // sp->AddPoint(isptrack,&p); isptrack++;
1472 // }
1473 }
1474 }
1475 }
1476 fTracker[3]->UnloadClusters();
1477 fLoader[3]->UnloadRecPoints();
1478}