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