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