git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding conventions
[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
c757bafd 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
596a855f 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
b26c3770 485//_____________________________________________________________________________
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
596a855f 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
1f46a9ae 613//_____________________________________________________________________________
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
2257f27e 673//_____________________________________________________________________________
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;
7b61cd9c 721 fLoader[iDet]->LoadRecPoints("read");
722 tree = fLoader[iDet]->TreeR();
b8cd5251 723 if (!tree) {
724 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 725 return kFALSE;
726 }
b8cd5251 727 fTracker[iDet]->LoadClusters(tree);
728 }
24f7a148 729
b8cd5251 730 // run tracking
731 if (fTracker[iDet]->PropagateBack(esd) != 0) {
732 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
733 return kFALSE;
734 }
735 if (fCheckPointLevel > 1) {
736 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
737 }
24f7a148 738
b8cd5251 739 // unload clusters
740 if (iDet > 2) { // all except ITS, TPC, TRD
741 fTracker[iDet]->UnloadClusters();
7b61cd9c 742 fLoader[iDet]->UnloadRecPoints();
b8cd5251 743 }
744 }
596a855f 745
b8cd5251 746 // pass 3: TRD + TPC + ITS refit inwards
747 for (Int_t iDet = 2; iDet >= 0; iDet--) {
748 if (!fTracker[iDet]) continue;
749 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 750
b8cd5251 751 // run tracking
752 if (fTracker[iDet]->RefitInward(esd) != 0) {
753 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
754 return kFALSE;
755 }
756 if (fCheckPointLevel > 1) {
757 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
758 }
596a855f 759
b8cd5251 760 // unload clusters
761 fTracker[iDet]->UnloadClusters();
762 fLoader[iDet]->UnloadRecPoints();
763 }
596a855f 764
5f8272e1 765 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
766 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 767
596a855f 768 return kTRUE;
769}
770
771//_____________________________________________________________________________
24f7a148 772Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 773{
774// fill the event summary data
775
030b532d 776 TStopwatch stopwatch;
777 stopwatch.Start();
815c2b38 778 AliInfo("filling ESD");
030b532d 779
596a855f 780 TString detStr = detectors;
b8cd5251 781 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
782 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
783 AliReconstructor* reconstructor = GetReconstructor(iDet);
784 if (!reconstructor) continue;
785
786 if (!ReadESD(esd, fgkDetectorName[iDet])) {
787 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 788 TTree* clustersTree = NULL;
789 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
790 fLoader[iDet]->LoadRecPoints("read");
791 clustersTree = fLoader[iDet]->TreeR();
792 if (!clustersTree) {
793 AliError(Form("Can't get the %s clusters tree",
794 fgkDetectorName[iDet]));
795 if (fStopOnError) return kFALSE;
796 }
797 }
798 if (fRawReader && !reconstructor->HasDigitConversion()) {
799 reconstructor->FillESD(fRawReader, clustersTree, esd);
800 } else {
801 TTree* digitsTree = NULL;
802 if (fLoader[iDet]) {
803 fLoader[iDet]->LoadDigits("read");
804 digitsTree = fLoader[iDet]->TreeD();
805 if (!digitsTree) {
806 AliError(Form("Can't get the %s digits tree",
807 fgkDetectorName[iDet]));
808 if (fStopOnError) return kFALSE;
809 }
810 }
811 reconstructor->FillESD(digitsTree, clustersTree, esd);
812 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
813 }
814 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
815 fLoader[iDet]->UnloadRecPoints();
816 }
817
b8cd5251 818 if (fRawReader) {
819 reconstructor->FillESD(fRunLoader, fRawReader, esd);
820 } else {
821 reconstructor->FillESD(fRunLoader, esd);
24f7a148 822 }
b8cd5251 823 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 824 }
825 }
826
827 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 828 AliError(Form("the following detectors were not found: %s",
829 detStr.Data()));
596a855f 830 if (fStopOnError) return kFALSE;
831 }
832
5f8272e1 833 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
834 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 835
596a855f 836 return kTRUE;
837}
838
839
840//_____________________________________________________________________________
841Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
842{
843// check whether detName is contained in detectors
844// if yes, it is removed from detectors
845
846 // check if all detectors are selected
847 if ((detectors.CompareTo("ALL") == 0) ||
848 detectors.BeginsWith("ALL ") ||
849 detectors.EndsWith(" ALL") ||
850 detectors.Contains(" ALL ")) {
851 detectors = "ALL";
852 return kTRUE;
853 }
854
855 // search for the given detector
856 Bool_t result = kFALSE;
857 if ((detectors.CompareTo(detName) == 0) ||
858 detectors.BeginsWith(detName+" ") ||
859 detectors.EndsWith(" "+detName) ||
860 detectors.Contains(" "+detName+" ")) {
861 detectors.ReplaceAll(detName, "");
862 result = kTRUE;
863 }
864
865 // clean up the detectors string
866 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
867 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
868 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
869
870 return result;
871}
e583c30d 872
f08fc9f5 873//_____________________________________________________________________________
874Bool_t AliReconstruction::InitRunLoader()
875{
876// get or create the run loader
877
878 if (gAlice) delete gAlice;
879 gAlice = NULL;
880
b26c3770 881 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
882 // load all base libraries to get the loader classes
883 TString libs = gSystem->GetLibraries();
884 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
885 TString detName = fgkDetectorName[iDet];
886 if (detName == "HLT") continue;
887 if (libs.Contains("lib" + detName + "base.so")) continue;
888 gSystem->Load("lib" + detName + "base.so");
889 }
f08fc9f5 890 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
891 if (!fRunLoader) {
892 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
893 CleanUp();
894 return kFALSE;
895 }
b26c3770 896 fRunLoader->CdGAFile();
897 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
898 if (fRunLoader->LoadgAlice() == 0) {
899 gAlice = fRunLoader->GetAliRun();
c84a5e9e 900 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
7ff9d95a 901 AliExternalTrackParam::SetFieldMap(gAlice->Field());
902 if(fUniformField)
903 AliExternalTrackParam::SetUniformFieldTracking();
904 else
905 AliExternalTrackParam::SetNonuniformFieldTracking();
b26c3770 906 }
f08fc9f5 907 }
908 if (!gAlice && !fRawReader) {
909 AliError(Form("no gAlice object found in file %s",
910 fGAliceFileName.Data()));
911 CleanUp();
912 return kFALSE;
913 }
914
915 } else { // galice.root does not exist
916 if (!fRawReader) {
917 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
918 CleanUp();
919 return kFALSE;
920 }
921 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
922 AliConfig::GetDefaultEventFolderName(),
923 "recreate");
924 if (!fRunLoader) {
925 AliError(Form("could not create run loader in file %s",
926 fGAliceFileName.Data()));
927 CleanUp();
928 return kFALSE;
929 }
930 fRunLoader->MakeTree("E");
931 Int_t iEvent = 0;
932 while (fRawReader->NextEvent()) {
933 fRunLoader->SetEventNumber(iEvent);
934 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
935 iEvent, iEvent);
936 fRunLoader->MakeTree("H");
937 fRunLoader->TreeE()->Fill();
938 iEvent++;
939 }
940 fRawReader->RewindEvents();
941 fRunLoader->WriteHeader("OVERWRITE");
942 fRunLoader->CdGAFile();
943 fRunLoader->Write(0, TObject::kOverwrite);
944// AliTracker::SetFieldMap(???);
945 }
946
947 return kTRUE;
948}
949
c757bafd 950//_____________________________________________________________________________
b8cd5251 951AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 952{
f08fc9f5 953// get the reconstructor object and the loader for a detector
c757bafd 954
b8cd5251 955 if (fReconstructor[iDet]) return fReconstructor[iDet];
956
957 // load the reconstructor object
958 TPluginManager* pluginManager = gROOT->GetPluginManager();
959 TString detName = fgkDetectorName[iDet];
960 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 961 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 962
963 if (detName == "HLT") {
964 if (!gROOT->GetClass("AliLevel3")) {
965 gSystem->Load("libAliL3Src.so");
966 gSystem->Load("libAliL3Misc.so");
967 gSystem->Load("libAliL3Hough.so");
968 gSystem->Load("libAliL3Comp.so");
969 }
970 }
971
972 AliReconstructor* reconstructor = NULL;
973 // first check if a plugin is defined for the reconstructor
974 TPluginHandler* pluginHandler =
975 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 976 // if not, add a plugin for it
977 if (!pluginHandler) {
b8cd5251 978 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 979 TString libs = gSystem->GetLibraries();
980 if (libs.Contains("lib" + detName + "base.so") ||
981 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 982 pluginManager->AddHandler("AliReconstructor", detName,
983 recName, detName + "rec", recName + "()");
984 } else {
985 pluginManager->AddHandler("AliReconstructor", detName,
986 recName, detName, recName + "()");
c757bafd 987 }
b8cd5251 988 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
989 }
990 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
991 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 992 }
b8cd5251 993 if (reconstructor) {
994 TObject* obj = fOptions.FindObject(detName.Data());
995 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 996 reconstructor->Init(fRunLoader);
b8cd5251 997 fReconstructor[iDet] = reconstructor;
998 }
999
f08fc9f5 1000 // get or create the loader
1001 if (detName != "HLT") {
1002 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1003 if (!fLoader[iDet]) {
1004 AliConfig::Instance()
1005 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1006 detName, detName);
1007 // first check if a plugin is defined for the loader
1008 TPluginHandler* pluginHandler =
1009 pluginManager->FindHandler("AliLoader", detName);
1010 // if not, add a plugin for it
1011 if (!pluginHandler) {
1012 TString loaderName = "Ali" + detName + "Loader";
1013 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1014 pluginManager->AddHandler("AliLoader", detName,
1015 loaderName, detName + "base",
1016 loaderName + "(const char*, TFolder*)");
1017 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1018 }
1019 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1020 fLoader[iDet] =
1021 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1022 fRunLoader->GetEventFolder());
1023 }
1024 if (!fLoader[iDet]) { // use default loader
1025 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1026 }
1027 if (!fLoader[iDet]) {
1028 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1029 if (fStopOnError) return NULL;
f08fc9f5 1030 } else {
1031 fRunLoader->AddLoader(fLoader[iDet]);
1032 fRunLoader->CdGAFile();
1033 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1034 fRunLoader->Write(0, TObject::kOverwrite);
1035 }
1036 }
1037 }
1038
b8cd5251 1039 return reconstructor;
c757bafd 1040}
1041
2257f27e 1042//_____________________________________________________________________________
1043Bool_t AliReconstruction::CreateVertexer()
1044{
1045// create the vertexer
1046
b8cd5251 1047 fVertexer = NULL;
1048 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1049 if (itsReconstructor) {
b8cd5251 1050 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1051 }
b8cd5251 1052 if (!fVertexer) {
815c2b38 1053 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1054 if (fStopOnError) return kFALSE;
1055 }
1056
1057 return kTRUE;
1058}
1059
24f7a148 1060//_____________________________________________________________________________
b8cd5251 1061Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1062{
f08fc9f5 1063// create the trackers
24f7a148 1064
b8cd5251 1065 TString detStr = detectors;
1066 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1067 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1068 AliReconstructor* reconstructor = GetReconstructor(iDet);
1069 if (!reconstructor) continue;
1070 TString detName = fgkDetectorName[iDet];
1f46a9ae 1071 if (detName == "HLT") {
1072 fRunHLTTracking = kTRUE;
1073 continue;
1074 }
f08fc9f5 1075
1076 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1077 if (!fTracker[iDet] && (iDet < 7)) {
1078 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1079 if (fStopOnError) return kFALSE;
1080 }
1081 }
1082
24f7a148 1083 return kTRUE;
1084}
1085
e583c30d 1086//_____________________________________________________________________________
b26c3770 1087void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1088{
1089// delete trackers and the run loader and close and delete the file
1090
b8cd5251 1091 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1092 delete fReconstructor[iDet];
1093 fReconstructor[iDet] = NULL;
1094 fLoader[iDet] = NULL;
1095 delete fTracker[iDet];
1096 fTracker[iDet] = NULL;
1097 }
1098 delete fVertexer;
1099 fVertexer = NULL;
e583c30d 1100
1101 delete fRunLoader;
1102 fRunLoader = NULL;
b649205a 1103 delete fRawReader;
1104 fRawReader = NULL;
e583c30d 1105
1106 if (file) {
1107 file->Close();
1108 delete file;
1109 }
b26c3770 1110
1111 if (fileOld) {
1112 fileOld->Close();
1113 delete fileOld;
1114 gSystem->Unlink("AliESDs.old.root");
1115 }
e583c30d 1116}
24f7a148 1117
1118
1119//_____________________________________________________________________________
1120Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1121{
1122// read the ESD event from a file
1123
1124 if (!esd) return kFALSE;
1125 char fileName[256];
1126 sprintf(fileName, "ESD_%d.%d_%s.root",
1127 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1128 if (gSystem->AccessPathName(fileName)) return kFALSE;
1129
f3a97c86 1130 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1131 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1132 TFile* file = TFile::Open(fileName);
1133 if (!file || !file->IsOpen()) {
815c2b38 1134 AliError(Form("opening %s failed", fileName));
24f7a148 1135 delete file;
1136 return kFALSE;
1137 }
1138
1139 gROOT->cd();
1140 delete esd;
1141 esd = (AliESD*) file->Get("ESD");
1142 file->Close();
1143 delete file;
1144 return kTRUE;
1145}
1146
1147//_____________________________________________________________________________
1148void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1149{
1150// write the ESD event to a file
1151
1152 if (!esd) return;
1153 char fileName[256];
1154 sprintf(fileName, "ESD_%d.%d_%s.root",
1155 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1156
815c2b38 1157 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1158 TFile* file = TFile::Open(fileName, "recreate");
1159 if (!file || !file->IsOpen()) {
815c2b38 1160 AliError(Form("opening %s failed", fileName));
24f7a148 1161 } else {
1162 esd->Write("ESD");
1163 file->Close();
1164 }
1165 delete file;
1166}
f3a97c86 1167
1168
1169
1170
1171//_____________________________________________________________________________
1172void AliReconstruction::CreateTag(TFile* file)
1173{
4302e20f 1174 Int_t ntrack;
1175 Int_t NProtons, NKaons, NPions, NMuons, NElectrons;
1176 Int_t Npos, Nneg, Nneutr;
1177 Int_t NK0s, Nneutrons, Npi0s, Ngamas;
1178 Int_t Nch1GeV, Nch3GeV, Nch10GeV;
1179 Int_t Nmu1GeV, Nmu3GeV, Nmu10GeV;
1180 Int_t Nel1GeV, Nel3GeV, Nel10GeV;
1181 Float_t MaxPt = .0, MeanPt = .0, TotalP = .0;
1182
f3a97c86 1183 AliRunTag *tag = new AliRunTag();
1184 AliDetectorTag *detTag = new AliDetectorTag();
1185 AliEventTag *evTag = new AliEventTag();
1186 TTree ttag("T","A Tree with event tags");
1187 TBranch * btag = ttag.Branch("AliTAG", "AliRunTag", &tag);
1188 btag->SetCompressionLevel(9);
1189
1190 AliInfo(Form("Creating the tags......."));
1191
1192 if (!file || !file->IsOpen()) {
1193 AliError(Form("opening failed"));
1194 delete file;
1195 return ;
1196 }
1197
1198 TTree *t = (TTree*) file->Get("esdTree");
1199 TBranch * b = t->GetBranch("ESD");
1200 AliESD *esd = 0;
1201 b->SetAddress(&esd);
1202
1203 tag->SetRunId(esd->GetRunNumber());
1204
1205 Int_t firstEvent = 0,lastEvent = 0;
1206 Int_t i_NumberOfEvents = b->GetEntries();
1207 for (Int_t i_EventNumber = 0; i_EventNumber < i_NumberOfEvents; i_EventNumber++)
1208 {
4302e20f 1209 ntrack = 0;
1210 Npos = 0;
1211 Nneg = 0;
1212 Nneutr =0;
1213 NK0s = 0;
1214 Nneutrons = 0;
1215 Npi0s = 0;
1216 Ngamas = 0;
1217 NProtons = 0;
1218 NKaons = 0;
1219 NPions = 0;
1220 NMuons = 0;
1221 NElectrons = 0;
1222 Nch1GeV = 0;
1223 Nch3GeV = 0;
1224 Nch10GeV = 0;
1225 Nmu1GeV = 0;
1226 Nmu3GeV = 0;
1227 Nmu10GeV = 0;
1228 Nel1GeV = 0;
1229 Nel3GeV = 0;
1230 Nel10GeV = 0;
1231 MaxPt = .0;
1232 MeanPt = .0;
1233 TotalP = .0;
1234
f3a97c86 1235 b->GetEntry(i_EventNumber);
f3a97c86 1236 const AliESDVertex * VertexIn = esd->GetVertex();
4302e20f 1237
1238 for (Int_t i_TrackNumber = 0; i_TrackNumber < esd->GetNumberOfTracks(); i_TrackNumber++)
1239 {
1240 AliESDtrack * ESDTrack = esd->GetTrack(i_TrackNumber);
1241 UInt_t status = ESDTrack->GetStatus();
1242
1243 //select only tracks with ITS refit
1244 if ((status&AliESDtrack::kITSrefit)==0) continue;
1245
1246 //select only tracks with TPC refit-->remove extremely high Pt tracks
1247 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1248
1249 //select only tracks with the "combined PID"
1250 if ((status&AliESDtrack::kESDpid)==0) continue;
1251 Double_t p[3];
1252 ESDTrack->GetPxPyPz(p);
1253 Double_t P = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1254 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1255 TotalP += P;
1256 MeanPt += fPt;
1257 if(fPt > MaxPt)
1258 MaxPt = fPt;
1259
1260 if(ESDTrack->GetSign() > 0)
1261 {
1262 Npos++;
1263 if(fPt > 1.0)
1264 Nch1GeV++;
1265 if(fPt > 3.0)
1266 Nch3GeV++;
1267 if(fPt > 10.0)
1268 Nch10GeV++;
1269 }
1270 if(ESDTrack->GetSign() < 0)
1271 {
1272 Nneg++;
1273 if(fPt > 1.0)
1274 Nch1GeV++;
1275 if(fPt > 3.0)
1276 Nch3GeV++;
1277 if(fPt > 10.0)
1278 Nch10GeV++;
1279 }
1280 if(ESDTrack->GetSign() == 0)
1281 Nneutr++;
1282
1283 //PID
1284 Double_t prob[10];
1285 ESDTrack->GetESDpid(prob);
1286
1287 //K0s
1288 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]))
1289 NK0s++;
1290 //neutrons
1291 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]))
1292 Nneutrons++;
1293 //pi0s
1294 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]))
1295 Npi0s++;
1296 //gamas
1297 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]))
1298 Ngamas++;
1299 //protons
1300 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]))
1301 NProtons++;
1302 //kaons
1303 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]))
1304 NKaons++;
1305 //kaons
1306 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]))
1307 NPions++;
1308 //muons
1309 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]))
1310 {
1311 NMuons++;
1312 if(fPt > 1.0)
1313 Nmu1GeV++;
1314 if(fPt > 3.0)
1315 Nmu3GeV++;
1316 if(fPt > 10.0)
1317 Nmu10GeV++;
1318 }
1319 //electrons
1320 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]))
1321 {
1322 NElectrons++;
1323 if(fPt > 1.0)
1324 Nel1GeV++;
1325 if(fPt > 3.0)
1326 Nel3GeV++;
1327 if(fPt > 10.0)
1328 Nel10GeV++;
1329 }
1330
1331
1332
1333 ntrack++;
1334 }//track loop
1335 // Fill the event tags
b085f7b9 1336 if(ntrack != 0)
1337 MeanPt = MeanPt/ntrack;
f3a97c86 1338
1339 evTag->SetEventId(i_EventNumber+1);
f3a97c86 1340 evTag->SetVertexX(VertexIn->GetXv());
1341 evTag->SetVertexY(VertexIn->GetYv());
1342 evTag->SetVertexZ(VertexIn->GetZv());
4302e20f 1343
1344 evTag->SetT0VertexZ(esd->GetT0zVertex());
1345
1346 evTag->SetTrigger(esd->GetTrigger());
1347
1348 evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
1349 evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
1350 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1351 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1352
1353
f3a97c86 1354 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
4302e20f 1355 evTag->SetNumOfPosTracks(Npos);
1356 evTag->SetNumOfNegTracks(Nneg);
1357 evTag->SetNumOfNeutrTracks(Nneutr);
1358
f3a97c86 1359 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1360 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
4302e20f 1361 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1362 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1363
1364 evTag->SetNumOfProtons(NProtons);
1365 evTag->SetNumOfKaons(NKaons);
1366 evTag->SetNumOfPions(NPions);
1367 evTag->SetNumOfMuons(NMuons);
1368 evTag->SetNumOfElectrons(NElectrons);
1369 evTag->SetNumOfPhotons(Ngamas);
1370 evTag->SetNumOfPi0s(Npi0s);
1371 evTag->SetNumOfNeutrons(Nneutrons);
1372 evTag->SetNumOfKaon0s(NK0s);
1373
1374 evTag->SetNumOfChargedAbove1GeV(Nch1GeV);
1375 evTag->SetNumOfChargedAbove3GeV(Nch3GeV);
1376 evTag->SetNumOfChargedAbove10GeV(Nch10GeV);
1377 evTag->SetNumOfMuonsAbove1GeV(Nmu1GeV);
1378 evTag->SetNumOfMuonsAbove3GeV(Nmu3GeV);
1379 evTag->SetNumOfMuonsAbove10GeV(Nmu10GeV);
1380 evTag->SetNumOfElectronsAbove1GeV(Nel1GeV);
1381 evTag->SetNumOfElectronsAbove3GeV(Nel3GeV);
1382 evTag->SetNumOfElectronsAbove10GeV(Nel10GeV);
f3a97c86 1383
1384 evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
1385 evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
4302e20f 1386
1387 evTag->SetTotalMomentum(TotalP);
1388 evTag->SetMeanPt(MeanPt);
1389 evTag->SetMaxPt(MaxPt);
1390
cb1645b7 1391 tag->AddEventTag(*evTag);
f3a97c86 1392 }
1393 lastEvent = i_NumberOfEvents;
1394
1395 ttag.Fill();
1396 tag->Clear();
1397
1398 char fileName[256];
1399 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1400 tag->GetRunId(),firstEvent,lastEvent );
1401 AliInfo(Form("writing tags to file %s", fileName));
1402 AliDebug(1, Form("writing tags to file %s", fileName));
1403
1404 TFile* ftag = TFile::Open(fileName, "recreate");
1405 ftag->cd();
1406 ttag.Write();
1407 ftag->Close();
1408 file->cd();
1409 delete tag;
1410 delete detTag;
1411 delete evTag;
1412}
1413