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