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