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