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