]> 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
Record changes.
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
CommitLineData
596a855f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// class for running the reconstruction //
21// //
22// Clusters and tracks are created for all detectors and all events by //
23// typing: //
24// //
25// AliReconstruction rec; //
26// rec.Run(); //
27// //
28// The Run method returns kTRUE in case of successful execution. //
c71de921 29// //
30// If the input to the reconstruction are not simulated digits but raw data, //
31// this can be specified by an argument of the Run method or by the method //
32// //
33// rec.SetInput("..."); //
34// //
35// The input formats and the corresponding argument are: //
36// - DDL raw data files: directory name, ends with "/" //
37// - raw data root file: root file name, extension ".root" //
38// - raw data DATE file: DATE file name, any other non-empty string //
39// - MC root files : empty string, default //
40// //
b26c3770 41// By default all events are reconstructed. The reconstruction can be //
42// limited to a range of events by giving the index of the first and the //
43// last event as an argument to the Run method or by calling //
44// //
45// rec.SetEventRange(..., ...); //
46// //
47// The index -1 (default) can be used for the last event to indicate no //
48// upper limit of the event range. //
49// //
973388c2 50// In case of raw-data reconstruction the user can modify the default //
51// number of events per digits/clusters/tracks file. In case the option //
52// is not used the number is set 1. In case the user provides 0, than //
53// the number of events is equal to the number of events inside the //
54// raw-data file (i.e. one digits/clusters/tracks file): //
55// //
56// rec.SetNumberOfEventsPerFile(...); //
57// //
58// //
596a855f 59// The name of the galice file can be changed from the default //
e583c30d 60// "galice.root" by passing it as argument to the AliReconstruction //
61// constructor or by //
596a855f 62// //
63// rec.SetGAliceFile("..."); //
64// //
59697224 65// The local reconstruction can be switched on or off for individual //
66// detectors by //
596a855f 67// //
59697224 68// rec.SetRunLocalReconstruction("..."); //
596a855f 69// //
70// The argument is a (case sensitive) string with the names of the //
71// detectors separated by a space. The special string "ALL" selects all //
72// available detectors. This is the default. //
73// //
c71de921 74// The reconstruction of the primary vertex position can be switched off by //
75// //
76// rec.SetRunVertexFinder(kFALSE); //
77// //
b8cd5251 78// The tracking and the creation of ESD tracks can be switched on for //
79// selected detectors by //
596a855f 80// //
b8cd5251 81// rec.SetRunTracking("..."); //
596a855f 82// //
c84a5e9e 83// Uniform/nonuniform field tracking switches (default: uniform field) //
84// //
1d99986f 85// rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
c84a5e9e 86// //
596a855f 87// The filling of additional ESD information can be steered by //
88// //
89// rec.SetFillESD("..."); //
90// //
b8cd5251 91// Again, for both methods the string specifies the list of detectors. //
92// The default is "ALL". //
93// //
94// The call of the shortcut method //
95// //
96// rec.SetRunReconstruction("..."); //
97// //
98// is equivalent to calling SetRunLocalReconstruction, SetRunTracking and //
99// SetFillESD with the same detector selecting string as argument. //
596a855f 100// //
c71de921 101// The reconstruction requires digits or raw data as input. For the creation //
102// of digits and raw data have a look at the class AliSimulation. //
596a855f 103// //
24f7a148 104// For debug purposes the method SetCheckPointLevel can be used. If the //
105// argument is greater than 0, files with ESD events will be written after //
106// selected steps of the reconstruction for each event: //
107// level 1: after tracking and after filling of ESD (final) //
108// level 2: in addition after each tracking step //
109// level 3: in addition after the filling of ESD for each detector //
110// If a final check point file exists for an event, this event will be //
111// skipped in the reconstruction. The tracking and the filling of ESD for //
112// a detector will be skipped as well, if the corresponding check point //
113// file exists. The ESD event will then be loaded from the file instead. //
114// //
596a855f 115///////////////////////////////////////////////////////////////////////////////
116
024a7e64 117#include <TArrayF.h>
118#include <TFile.h>
119#include <TSystem.h>
120#include <TROOT.h>
121#include <TPluginManager.h>
fd46e2d2 122#include <TStopwatch.h>
3103d196 123#include <TGeoManager.h>
2bdb9d38 124#include <TLorentzVector.h>
596a855f 125
126#include "AliReconstruction.h"
b8cd5251 127#include "AliReconstructor.h"
815c2b38 128#include "AliLog.h"
596a855f 129#include "AliRunLoader.h"
130#include "AliRun.h"
b649205a 131#include "AliRawReaderFile.h"
132#include "AliRawReaderDate.h"
133#include "AliRawReaderRoot.h"
001397cd 134#include "AliRawEventHeaderBase.h"
596a855f 135#include "AliESD.h"
1d99986f 136#include "AliESDfriend.h"
2257f27e 137#include "AliESDVertex.h"
32e449be 138#include "AliMultiplicity.h"
c84a5e9e 139#include "AliTracker.h"
2257f27e 140#include "AliVertexer.h"
c5e3e5d1 141#include "AliVertexerTracks.h"
5e4ff34d 142#include "AliV0vertexer.h"
143#include "AliCascadeVertexer.h"
596a855f 144#include "AliHeader.h"
145#include "AliGenEventHeader.h"
b26c3770 146#include "AliPID.h"
596a855f 147#include "AliESDpid.h"
ff8bb5ae 148#include "AliESDtrack.h"
f3a97c86 149
150#include "AliRunTag.h"
f3a97c86 151#include "AliDetectorTag.h"
152#include "AliEventTag.h"
153
98937d93 154#include "AliTrackPointArray.h"
b0314964 155#include "AliCDBManager.h"
6bae477a 156#include "AliCDBEntry.h"
157#include "AliAlignObj.h"
f3a97c86 158
b647652d 159#include "AliCentralTrigger.h"
160#include "AliCTPRawStream.h"
161
596a855f 162ClassImp(AliReconstruction)
163
164
c757bafd 165//_____________________________________________________________________________
b384f8a4 166const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
c757bafd 167
596a855f 168//_____________________________________________________________________________
024cf675 169AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 170 const char* name, const char* title) :
171 TNamed(name, title),
172
c84a5e9e 173 fUniformField(kTRUE),
2257f27e 174 fRunVertexFinder(kTRUE),
1f46a9ae 175 fRunHLTTracking(kFALSE),
e66fbafb 176 fRunMuonTracking(kFALSE),
1d99986f 177 fStopOnError(kFALSE),
178 fWriteAlignmentData(kFALSE),
179 fWriteESDfriend(kFALSE),
a7807689 180 fWriteAOD(kFALSE),
b647652d 181 fFillTriggerESD(kTRUE),
1d99986f 182
183 fRunLocalReconstruction("ALL"),
b8cd5251 184 fRunTracking("ALL"),
e583c30d 185 fFillESD("ALL"),
186 fGAliceFileName(gAliceFilename),
b649205a 187 fInput(""),
35042093 188 fEquipIdMap(""),
b26c3770 189 fFirstEvent(0),
190 fLastEvent(-1),
973388c2 191 fNumberOfEventsPerFile(1),
24f7a148 192 fCheckPointLevel(0),
b8cd5251 193 fOptions(),
6bae477a 194 fLoadAlignFromCDB(kTRUE),
195 fLoadAlignData("ALL"),
e583c30d 196
197 fRunLoader(NULL),
b649205a 198 fRawReader(NULL),
b8cd5251 199
98937d93 200 fVertexer(NULL),
9178838a 201 fDiamondProfile(NULL),
98937d93 202
6bae477a 203 fAlignObjArray(NULL),
ec92bee0 204 fCDBUri(cdbUri),
205 fSpecCDBUri()
596a855f 206{
207// create reconstruction object with default parameters
b8cd5251 208
209 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
210 fReconstructor[iDet] = NULL;
211 fLoader[iDet] = NULL;
212 fTracker[iDet] = NULL;
213 }
e47c4c2e 214 AliPID pid;
596a855f 215}
216
217//_____________________________________________________________________________
218AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 219 TNamed(rec),
220
c84a5e9e 221 fUniformField(rec.fUniformField),
2257f27e 222 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 223 fRunHLTTracking(rec.fRunHLTTracking),
e66fbafb 224 fRunMuonTracking(rec.fRunMuonTracking),
1d99986f 225 fStopOnError(rec.fStopOnError),
226 fWriteAlignmentData(rec.fWriteAlignmentData),
227 fWriteESDfriend(rec.fWriteESDfriend),
a7807689 228 fWriteAOD(rec.fWriteAOD),
b647652d 229 fFillTriggerESD(rec.fFillTriggerESD),
1d99986f 230
231 fRunLocalReconstruction(rec.fRunLocalReconstruction),
e583c30d 232 fRunTracking(rec.fRunTracking),
233 fFillESD(rec.fFillESD),
234 fGAliceFileName(rec.fGAliceFileName),
b649205a 235 fInput(rec.fInput),
35042093 236 fEquipIdMap(rec.fEquipIdMap),
b26c3770 237 fFirstEvent(rec.fFirstEvent),
238 fLastEvent(rec.fLastEvent),
973388c2 239 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
24f7a148 240 fCheckPointLevel(0),
b8cd5251 241 fOptions(),
6bae477a 242 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
243 fLoadAlignData(rec.fLoadAlignData),
e583c30d 244
245 fRunLoader(NULL),
b649205a 246 fRawReader(NULL),
b8cd5251 247
98937d93 248 fVertexer(NULL),
9178838a 249 fDiamondProfile(NULL),
98937d93 250
6bae477a 251 fAlignObjArray(rec.fAlignObjArray),
ec92bee0 252 fCDBUri(rec.fCDBUri),
253 fSpecCDBUri()
596a855f 254{
255// copy constructor
256
ec92bee0 257 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
efd2085e 258 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
259 }
b8cd5251 260 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
261 fReconstructor[iDet] = NULL;
262 fLoader[iDet] = NULL;
263 fTracker[iDet] = NULL;
264 }
ec92bee0 265 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
266 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
267 }
596a855f 268}
269
270//_____________________________________________________________________________
271AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
272{
273// assignment operator
274
275 this->~AliReconstruction();
276 new(this) AliReconstruction(rec);
277 return *this;
278}
279
280//_____________________________________________________________________________
281AliReconstruction::~AliReconstruction()
282{
283// clean up
284
e583c30d 285 CleanUp();
efd2085e 286 fOptions.Delete();
ec92bee0 287 fSpecCDBUri.Delete();
596a855f 288}
289
024cf675 290//_____________________________________________________________________________
291void AliReconstruction::InitCDBStorage()
292{
293// activate a default CDB storage
294// First check if we have any CDB storage set, because it is used
295// to retrieve the calibration and alignment constants
296
297 AliCDBManager* man = AliCDBManager::Instance();
ec92bee0 298 if (man->IsDefaultStorageSet())
024cf675 299 {
300 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 301 AliWarning("Default CDB storage has been already set !");
302 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
024cf675 303 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 304 fCDBUri = "";
305 }
306 else {
b8ec52f6 307 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
308 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
309 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 310 man->SetDefaultStorage(fCDBUri);
311 }
312
313 // Now activate the detector specific CDB storage locations
c3a7b59a 314 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
315 TObject* obj = fSpecCDBUri[i];
316 if (!obj) continue;
b8ec52f6 317 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
318 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
319 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
c3a7b59a 320 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
ec92bee0 321 }
727d0766 322 man->Print();
024cf675 323}
324
325//_____________________________________________________________________________
326void AliReconstruction::SetDefaultStorage(const char* uri) {
ec92bee0 327// Store the desired default CDB storage location
328// Activate it later within the Run() method
024cf675 329
ec92bee0 330 fCDBUri = uri;
024cf675 331
332}
333
334//_____________________________________________________________________________
c3a7b59a 335void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
ec92bee0 336// Store a detector-specific CDB storage location
337// Activate it later within the Run() method
024cf675 338
c3a7b59a 339 AliCDBPath aPath(calibType);
340 if(!aPath.IsValid()){
341 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
342 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
343 if(!strcmp(calibType, fgkDetectorName[iDet])) {
344 aPath.SetPath(Form("%s/*", calibType));
345 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
346 break;
347 }
348 }
349 if(!aPath.IsValid()){
350 AliError(Form("Not a valid path or detector: %s", calibType));
351 return;
352 }
353 }
354
355 // check that calibType refers to a "valid" detector name
356 Bool_t isDetector = kFALSE;
357 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
358 TString detName = fgkDetectorName[iDet];
359 if(aPath.GetLevel0() == detName) {
360 isDetector = kTRUE;
361 break;
362 }
363 }
364
365 if(!isDetector) {
366 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
367 return;
368 }
369
370 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
ec92bee0 371 if (obj) fSpecCDBUri.Remove(obj);
c3a7b59a 372 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
024cf675 373
374}
375
6bae477a 376//_____________________________________________________________________________
377Bool_t AliReconstruction::SetRunNumber()
378{
379 // The method is called in Run() in order
380 // to set a correct run number.
381 // In case of raw data reconstruction the
382 // run number is taken from the raw data header
383
384 if(AliCDBManager::Instance()->GetRun() < 0) {
385 if (!fRunLoader) {
386 AliError("No run loader is found !");
387 return kFALSE;
388 }
389 // read run number from gAlice
ec92bee0 390 if(fRunLoader->GetAliRun())
391 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
392 else {
393 if(fRawReader) {
394 if(fRawReader->NextEvent()) {
395 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
396 fRawReader->RewindEvents();
397 }
398 else {
399 AliError("No raw-data events found !");
400 return kFALSE;
401 }
402 }
403 else {
404 AliError("Neither gAlice nor RawReader objects are found !");
405 return kFALSE;
406 }
407 }
408 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
6bae477a 409 }
410 return kTRUE;
411}
412
413//_____________________________________________________________________________
414Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
415{
416 // Read collection of alignment objects (AliAlignObj derived) saved
417 // in the TClonesArray ClArrayName and apply them to the geometry
418 // manager singleton.
419 //
420 alObjArray->Sort();
421 Int_t nvols = alObjArray->GetEntriesFast();
422
dc0984f8 423 Bool_t flag = kTRUE;
424
6bae477a 425 for(Int_t j=0; j<nvols; j++)
426 {
427 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
dc0984f8 428 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
6bae477a 429 }
430
431 if (AliDebugLevelClass() >= 1) {
e0625db4 432 gGeoManager->GetTopNode()->CheckOverlaps(20);
6bae477a 433 TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
434 if(ovexlist->GetEntriesFast()){
435 AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
436 }
437 }
438
dc0984f8 439 return flag;
6bae477a 440
441}
442
443//_____________________________________________________________________________
444Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
445{
446 // Fills array of single detector's alignable objects from CDB
447
448 AliDebug(2, Form("Loading alignment data for detector: %s",detName));
449
450 AliCDBEntry *entry;
451
452 AliCDBPath path(detName,"Align","Data");
453
454 entry=AliCDBManager::Instance()->Get(path.GetPath());
455 if(!entry){
456 AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
457 return kFALSE;
458 }
459 entry->SetOwner(1);
460 TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
461 alignArray->SetOwner(0);
462 AliDebug(2,Form("Found %d alignment objects for %s",
463 alignArray->GetEntries(),detName));
464
465 AliAlignObj *alignObj=0;
466 TIter iter(alignArray);
467
468 // loop over align objects in detector
469 while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
470 fAlignObjArray->Add(alignObj);
471 }
472 // delete entry --- Don't delete, it is cached!
473
474 AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
475 return kTRUE;
476
477}
478
479//_____________________________________________________________________________
480Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
481{
482 // Read the alignment objects from CDB.
483 // Each detector is supposed to have the
484 // alignment objects in DET/Align/Data CDB path.
485 // All the detector objects are then collected,
486 // sorted by geometry level (starting from ALIC) and
487 // then applied to the TGeo geometry.
488 // Finally an overlaps check is performed.
489
490 // Load alignment data from CDB and fill fAlignObjArray
491 if(fLoadAlignFromCDB){
492 if(!fAlignObjArray) fAlignObjArray = new TObjArray();
493
494 //fAlignObjArray->RemoveAll();
495 fAlignObjArray->Clear();
496 fAlignObjArray->SetOwner(0);
497
498 TString detStr = detectors;
499 TString dataNotLoaded="";
500 TString dataLoaded="";
501
502 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
503 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
504 if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
505 dataNotLoaded += fgkDetectorName[iDet];
506 dataNotLoaded += " ";
507 } else {
508 dataLoaded += fgkDetectorName[iDet];
509 dataLoaded += " ";
510 }
511 } // end loop over detectors
512
513 if ((detStr.CompareTo("ALL") == 0)) detStr = "";
514 dataNotLoaded += detStr;
7250c343 515 if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
6bae477a 516 dataLoaded.Data()));
7250c343 517 if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
6bae477a 518 dataNotLoaded.Data()));
519 } // fLoadAlignFromCDB flag
520
521 // Check if the array with alignment objects was
522 // provided by the user. If yes, apply the objects
523 // to the present TGeo geometry
524 if (fAlignObjArray) {
525 if (gGeoManager && gGeoManager->IsClosed()) {
526 if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
dc0984f8 527 AliError("The misalignment of one or more volumes failed!"
528 "Compare the list of simulated detectors and the list of detector alignment data!");
6bae477a 529 return kFALSE;
530 }
531 }
532 else {
533 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
534 return kFALSE;
535 }
536 }
537
8e245d15 538 delete fAlignObjArray; fAlignObjArray=0;
6bae477a 539 return kTRUE;
540}
596a855f 541
542//_____________________________________________________________________________
543void AliReconstruction::SetGAliceFile(const char* fileName)
544{
545// set the name of the galice file
546
547 fGAliceFileName = fileName;
548}
549
efd2085e 550//_____________________________________________________________________________
551void AliReconstruction::SetOption(const char* detector, const char* option)
552{
553// set options for the reconstruction of a detector
554
555 TObject* obj = fOptions.FindObject(detector);
556 if (obj) fOptions.Remove(obj);
557 fOptions.Add(new TNamed(detector, option));
558}
559
596a855f 560
561//_____________________________________________________________________________
4a33489c 562Bool_t AliReconstruction::Run(const char* input)
596a855f 563{
564// run the reconstruction
565
b649205a 566 // set the input
567 if (!input) input = fInput.Data();
568 TString fileName(input);
569 if (fileName.EndsWith("/")) {
570 fRawReader = new AliRawReaderFile(fileName);
571 } else if (fileName.EndsWith(".root")) {
572 fRawReader = new AliRawReaderRoot(fileName);
573 } else if (!fileName.IsNull()) {
574 fRawReader = new AliRawReaderDate(fileName);
575 fRawReader->SelectEvents(7);
576 }
35042093 577 if (!fEquipIdMap.IsNull() && fRawReader)
578 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
579
b649205a 580
f08fc9f5 581 // get the run loader
582 if (!InitRunLoader()) return kFALSE;
596a855f 583
ec92bee0 584 // Initialize the CDB storage
585 InitCDBStorage();
586
6bae477a 587 // Set run number in CDBManager (if it is not already set by the user)
588 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
589
590 // Import ideal TGeo geometry and apply misalignment
591 if (!gGeoManager) {
592 TString geom(gSystem->DirName(fGAliceFileName));
593 geom += "/geometry.root";
594 TGeoManager::Import(geom.Data());
595 if (!gGeoManager) if (fStopOnError) return kFALSE;
596 }
8e245d15 597
598 AliCDBManager* man = AliCDBManager::Instance();
6bae477a 599 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
600
596a855f 601 // local reconstruction
59697224 602 if (!fRunLocalReconstruction.IsNull()) {
603 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 604 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 605 }
606 }
b26c3770 607// if (!fRunVertexFinder && fRunTracking.IsNull() &&
608// fFillESD.IsNull()) return kTRUE;
2257f27e 609
610 // get vertexer
611 if (fRunVertexFinder && !CreateVertexer()) {
612 if (fStopOnError) {
613 CleanUp();
614 return kFALSE;
615 }
616 }
596a855f 617
f08fc9f5 618 // get trackers
b8cd5251 619 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 620 if (fStopOnError) {
621 CleanUp();
622 return kFALSE;
623 }
596a855f 624 }
24f7a148 625
9db3a215 626
627 TStopwatch stopwatch;
628 stopwatch.Start();
629
b26c3770 630 // get the possibly already existing ESD file and tree
1f46a9ae 631 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
b26c3770 632 TFile* fileOld = NULL;
1f46a9ae 633 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 634 if (!gSystem->AccessPathName("AliESDs.root")){
635 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
636 fileOld = TFile::Open("AliESDs.old.root");
637 if (fileOld && fileOld->IsOpen()) {
638 treeOld = (TTree*) fileOld->Get("esdTree");
639 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
1f46a9ae 640 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
641 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
b26c3770 642 }
643 }
644
36711aa4 645 // create the ESD output file and tree
596a855f 646 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
647 if (!file->IsOpen()) {
815c2b38 648 AliError("opening AliESDs.root failed");
b26c3770 649 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 650 }
36711aa4 651 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
652 tree->Branch("ESD", "AliESD", &esd);
1f46a9ae 653 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
654 hlttree->Branch("ESD", "AliESD", &hltesd);
655 delete esd; delete hltesd;
656 esd = NULL; hltesd = NULL;
1d99986f 657
500d54ab 658 // create the branch with ESD additions
659 AliESDfriend *esdf=0;
1d99986f 660 if (fWriteESDfriend) {
500d54ab 661 TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
662 br->SetFile("AliESDfriends.root");
1d99986f 663 }
664
c5e3e5d1 665 AliVertexerTracks tVertexer;
9178838a 666 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
c5e3e5d1 667
596a855f 668 // loop over events
b649205a 669 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 670
596a855f 671 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 672 if (fRawReader) fRawReader->NextEvent();
4a33489c 673 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
b26c3770 674 // copy old ESD to the new one
675 if (treeOld) {
676 treeOld->SetBranchAddress("ESD", &esd);
677 treeOld->GetEntry(iEvent);
678 }
679 tree->Fill();
1f46a9ae 680 if (hlttreeOld) {
681 hlttreeOld->SetBranchAddress("ESD", &hltesd);
682 hlttreeOld->GetEntry(iEvent);
683 }
684 hlttree->Fill();
b26c3770 685 continue;
686 }
687
815c2b38 688 AliInfo(Form("processing event %d", iEvent));
596a855f 689 fRunLoader->GetEvent(iEvent);
24f7a148 690
bb0901a4 691 char aFileName[256];
692 sprintf(aFileName, "ESD_%d.%d_final.root",
f08fc9f5 693 fRunLoader->GetHeader()->GetRun(),
694 fRunLoader->GetHeader()->GetEventNrInRun());
bb0901a4 695 if (!gSystem->AccessPathName(aFileName)) continue;
24f7a148 696
b26c3770 697 // local reconstruction
698 if (!fRunLocalReconstruction.IsNull()) {
699 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
700 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
701 }
702 }
703
1f46a9ae 704 esd = new AliESD; hltesd = new AliESD;
f08fc9f5 705 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 706 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
31fd97b2 707 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
708 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
d6ee376f 709
710 // Set magnetic field from the tracker
711 esd->SetMagneticField(AliTracker::GetBz());
712 hltesd->SetMagneticField(AliTracker::GetBz());
596a855f 713
2e3550da 714 // Fill raw-data error log into the ESD
715 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
716
2257f27e 717 // vertex finder
718 if (fRunVertexFinder) {
719 if (!ReadESD(esd, "vertex")) {
720 if (!RunVertexFinder(esd)) {
b26c3770 721 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 722 }
723 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
724 }
725 }
726
1f46a9ae 727 // HLT tracking
728 if (!fRunTracking.IsNull()) {
729 if (fRunHLTTracking) {
730 hltesd->SetVertex(esd->GetVertex());
731 if (!RunHLTTracking(hltesd)) {
732 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
733 }
734 }
735 }
736
e66fbafb 737 // Muon tracking
b8cd5251 738 if (!fRunTracking.IsNull()) {
e66fbafb 739 if (fRunMuonTracking) {
740 if (!RunMuonTracking()) {
b26c3770 741 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 742 }
596a855f 743 }
744 }
745
e66fbafb 746 // barrel tracking
747 if (!fRunTracking.IsNull()) {
21c573b7 748 if (!ReadESD(esd, "tracking")) {
749 if (!RunTracking(esd)) {
750 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
e66fbafb 751 }
21c573b7 752 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
e66fbafb 753 }
754 }
21c573b7 755
596a855f 756 // fill ESD
757 if (!fFillESD.IsNull()) {
758 if (!FillESD(esd, fFillESD)) {
b26c3770 759 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 760 }
761 }
001397cd 762 // fill Event header information from the RawEventHeader
763 if (fRawReader){FillRawEventHeaderESD(esd);}
596a855f 764
765 // combined PID
766 AliESDpid::MakePID(esd);
24f7a148 767 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 768
b647652d 769 if (fFillTriggerESD) {
770 if (!ReadESD(esd, "trigger")) {
771 if (!FillTriggerESD(esd)) {
772 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
773 }
774 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
775 }
776 }
777
a03cd2e0 778 //Try to improve the reconstructed primary vertex position using the tracks
779 AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
780 if (pvtx)
781 if (pvtx->GetStatus()) {
782 // Store the improved primary vertex
783 esd->SetPrimaryVertex(pvtx);
784 // Propagate the tracks to the DCA to the improved primary vertex
785 Double_t somethingbig = 777.;
786 Double_t bz = esd->GetMagneticField();
787 Int_t nt=esd->GetNumberOfTracks();
788 while (nt--) {
789 AliESDtrack *t = esd->GetTrack(nt);
790 t->RelateToVertex(pvtx, bz, somethingbig);
791 }
792 }
c5e3e5d1 793
5e4ff34d 794 {
795 // V0 finding
796 AliV0vertexer vtxer;
797 vtxer.Tracks2V0vertices(esd);
798
799 // Cascade finding
800 AliCascadeVertexer cvtxer;
801 cvtxer.V0sTracks2CascadeVertices(esd);
802 }
803
596a855f 804 // write ESD
1d99986f 805 if (fWriteESDfriend) {
d75007f6 806 esdf=new AliESDfriend();
807 esd->GetESDfriend(esdf);
1d99986f 808 }
500d54ab 809 tree->Fill();
810
811 // write HLT ESD
812 hlttree->Fill();
1d99986f 813
f3a97c86 814 if (fCheckPointLevel > 0) WriteESD(esd, "final");
815
500d54ab 816 delete esd; delete esdf; delete hltesd;
817 esd = NULL; esdf=NULL; hltesd = NULL;
596a855f 818 }
819
9db3a215 820 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
821 stopwatch.RealTime(),stopwatch.CpuTime()));
822
36711aa4 823 file->cd();
a9c0e6db 824 if (fWriteESDfriend)
825 tree->SetBranchStatus("ESDfriend*",0);
36711aa4 826 tree->Write();
1f46a9ae 827 hlttree->Write();
f3a97c86 828
a7807689 829 if (fWriteAOD) {
830 CreateAOD(file);
831 }
832
f3a97c86 833 // Create tags for the events in the ESD tree (the ESD tree is always present)
834 // In case of empty events the tags will contain dummy values
835 CreateTag(file);
b26c3770 836 CleanUp(file, fileOld);
596a855f 837
838 return kTRUE;
839}
840
841
842//_____________________________________________________________________________
59697224 843Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 844{
59697224 845// run the local reconstruction
596a855f 846
030b532d 847 TStopwatch stopwatch;
848 stopwatch.Start();
849
8e245d15 850 AliCDBManager* man = AliCDBManager::Instance();
851 Bool_t origCache = man->GetCacheFlag();
852
596a855f 853 TString detStr = detectors;
b8cd5251 854 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
855 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
856 AliReconstructor* reconstructor = GetReconstructor(iDet);
857 if (!reconstructor) continue;
b26c3770 858 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 859
860 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
861 TStopwatch stopwatchDet;
862 stopwatchDet.Start();
8e245d15 863
864 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
865
866 man->SetCacheFlag(kTRUE);
867 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
868 man->GetAll(calibPath); // entries are cached!
869
b8cd5251 870 if (fRawReader) {
871 fRawReader->RewindEvents();
872 reconstructor->Reconstruct(fRunLoader, fRawReader);
873 } else {
874 reconstructor->Reconstruct(fRunLoader);
596a855f 875 }
5f8272e1 876 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
877 fgkDetectorName[iDet],
878 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
8e245d15 879
880 // unload calibration data
881 man->ClearCache();
596a855f 882 }
883
8e245d15 884 man->SetCacheFlag(origCache);
885
596a855f 886 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 887 AliError(Form("the following detectors were not found: %s",
888 detStr.Data()));
596a855f 889 if (fStopOnError) return kFALSE;
890 }
891
5f8272e1 892 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
893 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 894
596a855f 895 return kTRUE;
896}
897
b26c3770 898//_____________________________________________________________________________
899Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
900{
901// run the local reconstruction
902
903 TStopwatch stopwatch;
904 stopwatch.Start();
905
906 TString detStr = detectors;
907 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
908 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
909 AliReconstructor* reconstructor = GetReconstructor(iDet);
910 if (!reconstructor) continue;
911 AliLoader* loader = fLoader[iDet];
912
913 // conversion of digits
914 if (fRawReader && reconstructor->HasDigitConversion()) {
915 AliInfo(Form("converting raw data digits into root objects for %s",
916 fgkDetectorName[iDet]));
917 TStopwatch stopwatchDet;
918 stopwatchDet.Start();
919 loader->LoadDigits("update");
920 loader->CleanDigits();
921 loader->MakeDigitsContainer();
922 TTree* digitsTree = loader->TreeD();
923 reconstructor->ConvertDigits(fRawReader, digitsTree);
924 loader->WriteDigits("OVERWRITE");
925 loader->UnloadDigits();
5f8272e1 926 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
927 fgkDetectorName[iDet],
928 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 929 }
930
931 // local reconstruction
932 if (!reconstructor->HasLocalReconstruction()) continue;
933 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
934 TStopwatch stopwatchDet;
935 stopwatchDet.Start();
936 loader->LoadRecPoints("update");
937 loader->CleanRecPoints();
938 loader->MakeRecPointsContainer();
939 TTree* clustersTree = loader->TreeR();
940 if (fRawReader && !reconstructor->HasDigitConversion()) {
941 reconstructor->Reconstruct(fRawReader, clustersTree);
942 } else {
943 loader->LoadDigits("read");
944 TTree* digitsTree = loader->TreeD();
945 if (!digitsTree) {
946 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
947 if (fStopOnError) return kFALSE;
948 } else {
949 reconstructor->Reconstruct(digitsTree, clustersTree);
950 }
951 loader->UnloadDigits();
952 }
953 loader->WriteRecPoints("OVERWRITE");
954 loader->UnloadRecPoints();
5f8272e1 955 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
956 fgkDetectorName[iDet],
957 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 958 }
959
960 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
961 AliError(Form("the following detectors were not found: %s",
962 detStr.Data()));
963 if (fStopOnError) return kFALSE;
964 }
5f8272e1 965
966 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
967 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 968
969 return kTRUE;
970}
971
596a855f 972//_____________________________________________________________________________
2257f27e 973Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 974{
975// run the barrel tracking
976
030b532d 977 TStopwatch stopwatch;
978 stopwatch.Start();
979
2257f27e 980 AliESDVertex* vertex = NULL;
981 Double_t vtxPos[3] = {0, 0, 0};
982 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
983 TArrayF mcVertex(3);
a6b0b91b 984 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
985 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
986 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
987 }
2257f27e 988
b8cd5251 989 if (fVertexer) {
815c2b38 990 AliInfo("running the ITS vertex finder");
b26c3770 991 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 992 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 993 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 994 if(!vertex){
815c2b38 995 AliWarning("Vertex not found");
c710f220 996 vertex = new AliESDVertex();
d1a50cb5 997 vertex->SetName("default");
2257f27e 998 }
999 else {
1000 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 1001 vertex->SetName("reconstructed");
2257f27e 1002 }
1003
1004 } else {
815c2b38 1005 AliInfo("getting the primary vertex from MC");
2257f27e 1006 vertex = new AliESDVertex(vtxPos, vtxErr);
1007 }
1008
1009 if (vertex) {
1010 vertex->GetXYZ(vtxPos);
1011 vertex->GetSigmaXYZ(vtxErr);
1012 } else {
815c2b38 1013 AliWarning("no vertex reconstructed");
2257f27e 1014 vertex = new AliESDVertex(vtxPos, vtxErr);
1015 }
1016 esd->SetVertex(vertex);
32e449be 1017 // if SPD multiplicity has been determined, it is stored in the ESD
1018 AliMultiplicity *mult= fVertexer->GetMultiplicity();
1019 if(mult)esd->SetMultiplicity(mult);
1020
b8cd5251 1021 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1022 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1023 }
2257f27e 1024 delete vertex;
1025
5f8272e1 1026 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1027 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 1028
1029 return kTRUE;
1030}
1031
1f46a9ae 1032//_____________________________________________________________________________
1033Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1034{
1035// run the HLT barrel tracking
1036
1037 TStopwatch stopwatch;
1038 stopwatch.Start();
1039
1040 if (!fRunLoader) {
1041 AliError("Missing runLoader!");
1042 return kFALSE;
1043 }
1044
1045 AliInfo("running HLT tracking");
1046
1047 // Get a pointer to the HLT reconstructor
1048 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1049 if (!reconstructor) return kFALSE;
1050
1051 // TPC + ITS
1052 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1053 TString detName = fgkDetectorName[iDet];
1054 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1055 reconstructor->SetOption(detName.Data());
1056 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1057 if (!tracker) {
1058 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1059 if (fStopOnError) return kFALSE;
9dcc06e1 1060 continue;
1f46a9ae 1061 }
1062 Double_t vtxPos[3];
1063 Double_t vtxErr[3]={0.005,0.005,0.010};
1064 const AliESDVertex *vertex = esd->GetVertex();
1065 vertex->GetXYZ(vtxPos);
1066 tracker->SetVertex(vtxPos,vtxErr);
1067 if(iDet != 1) {
1068 fLoader[iDet]->LoadRecPoints("read");
1069 TTree* tree = fLoader[iDet]->TreeR();
1070 if (!tree) {
1071 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1072 return kFALSE;
1073 }
1074 tracker->LoadClusters(tree);
1075 }
1076 if (tracker->Clusters2Tracks(esd) != 0) {
1077 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1078 return kFALSE;
1079 }
1080 if(iDet != 1) {
1081 tracker->UnloadClusters();
1082 }
1083 delete tracker;
1084 }
1085
5f8272e1 1086 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1087 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 1088
1089 return kTRUE;
1090}
1091
e66fbafb 1092//_____________________________________________________________________________
1093Bool_t AliReconstruction::RunMuonTracking()
1094{
1095// run the muon spectrometer tracking
1096
1097 TStopwatch stopwatch;
1098 stopwatch.Start();
1099
1100 if (!fRunLoader) {
1101 AliError("Missing runLoader!");
1102 return kFALSE;
1103 }
1104 Int_t iDet = 7; // for MUON
1105
1106 AliInfo("is running...");
1107
1108 // Get a pointer to the MUON reconstructor
1109 AliReconstructor *reconstructor = GetReconstructor(iDet);
1110 if (!reconstructor) return kFALSE;
1111
1112
1113 TString detName = fgkDetectorName[iDet];
1114 AliDebug(1, Form("%s tracking", detName.Data()));
1115 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1116 if (!tracker) {
1117 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1118 return kFALSE;
1119 }
1120
1121 // create Tracks
1122 fLoader[iDet]->LoadTracks("update");
1123 fLoader[iDet]->CleanTracks();
1124 fLoader[iDet]->MakeTracksContainer();
1125
1126 // read RecPoints
1127 fLoader[iDet]->LoadRecPoints("read");
1128
1129 if (!tracker->Clusters2Tracks(0x0)) {
1130 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1131 return kFALSE;
1132 }
1133 fLoader[iDet]->UnloadRecPoints();
1134
1135 fLoader[iDet]->WriteTracks("OVERWRITE");
1136 fLoader[iDet]->UnloadTracks();
1137
1138 delete tracker;
1139
1140
1141 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1142 stopwatch.RealTime(),stopwatch.CpuTime()));
1143
1144 return kTRUE;
1145}
1146
1147
2257f27e 1148//_____________________________________________________________________________
1149Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1150{
1151// run the barrel tracking
1152
1153 TStopwatch stopwatch;
1154 stopwatch.Start();
24f7a148 1155
815c2b38 1156 AliInfo("running tracking");
596a855f 1157
b8cd5251 1158 // pass 1: TPC + ITS inwards
1159 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1160 if (!fTracker[iDet]) continue;
1161 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 1162
b8cd5251 1163 // load clusters
1164 fLoader[iDet]->LoadRecPoints("read");
1165 TTree* tree = fLoader[iDet]->TreeR();
1166 if (!tree) {
1167 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1168 return kFALSE;
1169 }
b8cd5251 1170 fTracker[iDet]->LoadClusters(tree);
1171
1172 // run tracking
1173 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1174 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 1175 return kFALSE;
1176 }
b8cd5251 1177 if (fCheckPointLevel > 1) {
1178 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1179 }
878e1fe1 1180 // preliminary PID in TPC needed by the ITS tracker
1181 if (iDet == 1) {
1182 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 1183 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 1184 AliESDpid::MakePID(esd);
1185 }
b8cd5251 1186 }
596a855f 1187
b8cd5251 1188 // pass 2: ALL backwards
1189 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1190 if (!fTracker[iDet]) continue;
1191 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1192
1193 // load clusters
1194 if (iDet > 1) { // all except ITS, TPC
1195 TTree* tree = NULL;
7b61cd9c 1196 fLoader[iDet]->LoadRecPoints("read");
1197 tree = fLoader[iDet]->TreeR();
b8cd5251 1198 if (!tree) {
1199 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1200 return kFALSE;
1201 }
b8cd5251 1202 fTracker[iDet]->LoadClusters(tree);
1203 }
24f7a148 1204
b8cd5251 1205 // run tracking
1206 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1207 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1208 return kFALSE;
1209 }
1210 if (fCheckPointLevel > 1) {
1211 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1212 }
24f7a148 1213
b8cd5251 1214 // unload clusters
1215 if (iDet > 2) { // all except ITS, TPC, TRD
1216 fTracker[iDet]->UnloadClusters();
7b61cd9c 1217 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1218 }
8f37df88 1219 // updated PID in TPC needed by the ITS tracker -MI
1220 if (iDet == 1) {
1221 GetReconstructor(1)->FillESD(fRunLoader, esd);
1222 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1223 AliESDpid::MakePID(esd);
1224 }
b8cd5251 1225 }
596a855f 1226
98937d93 1227 // write space-points to the ESD in case alignment data output
1228 // is switched on
1229 if (fWriteAlignmentData)
1230 WriteAlignmentData(esd);
1231
b8cd5251 1232 // pass 3: TRD + TPC + ITS refit inwards
1233 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1234 if (!fTracker[iDet]) continue;
1235 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 1236
b8cd5251 1237 // run tracking
1238 if (fTracker[iDet]->RefitInward(esd) != 0) {
1239 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1240 return kFALSE;
1241 }
1242 if (fCheckPointLevel > 1) {
1243 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1244 }
596a855f 1245
b8cd5251 1246 // unload clusters
1247 fTracker[iDet]->UnloadClusters();
1248 fLoader[iDet]->UnloadRecPoints();
1249 }
ff8bb5ae 1250 //
1251 // Propagate track to the vertex - if not done by ITS
1252 //
1253 Int_t ntracks = esd->GetNumberOfTracks();
1254 for (Int_t itrack=0; itrack<ntracks; itrack++){
1255 const Double_t kRadius = 3; // beam pipe radius
1256 const Double_t kMaxStep = 5; // max step
1257 const Double_t kMaxD = 123456; // max distance to prim vertex
1258 Double_t fieldZ = AliTracker::GetBz(); //
1259 AliESDtrack * track = esd->GetTrack(itrack);
1260 if (!track) continue;
1261 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
a28195f7 1262 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
ff8bb5ae 1263 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1264 }
1265
5f8272e1 1266 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1267 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1268
596a855f 1269 return kTRUE;
1270}
1271
1272//_____________________________________________________________________________
24f7a148 1273Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1274{
1275// fill the event summary data
1276
030b532d 1277 TStopwatch stopwatch;
1278 stopwatch.Start();
815c2b38 1279 AliInfo("filling ESD");
030b532d 1280
596a855f 1281 TString detStr = detectors;
b8cd5251 1282 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1283 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1284 AliReconstructor* reconstructor = GetReconstructor(iDet);
1285 if (!reconstructor) continue;
1286
1287 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1288 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 1289 TTree* clustersTree = NULL;
1290 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1291 fLoader[iDet]->LoadRecPoints("read");
1292 clustersTree = fLoader[iDet]->TreeR();
1293 if (!clustersTree) {
1294 AliError(Form("Can't get the %s clusters tree",
1295 fgkDetectorName[iDet]));
1296 if (fStopOnError) return kFALSE;
1297 }
1298 }
1299 if (fRawReader && !reconstructor->HasDigitConversion()) {
1300 reconstructor->FillESD(fRawReader, clustersTree, esd);
1301 } else {
1302 TTree* digitsTree = NULL;
1303 if (fLoader[iDet]) {
1304 fLoader[iDet]->LoadDigits("read");
1305 digitsTree = fLoader[iDet]->TreeD();
1306 if (!digitsTree) {
1307 AliError(Form("Can't get the %s digits tree",
1308 fgkDetectorName[iDet]));
1309 if (fStopOnError) return kFALSE;
1310 }
1311 }
1312 reconstructor->FillESD(digitsTree, clustersTree, esd);
1313 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1314 }
1315 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1316 fLoader[iDet]->UnloadRecPoints();
1317 }
1318
b8cd5251 1319 if (fRawReader) {
1320 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1321 } else {
1322 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1323 }
b8cd5251 1324 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1325 }
1326 }
1327
1328 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1329 AliError(Form("the following detectors were not found: %s",
1330 detStr.Data()));
596a855f 1331 if (fStopOnError) return kFALSE;
1332 }
1333
5f8272e1 1334 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1335 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1336
596a855f 1337 return kTRUE;
1338}
1339
b647652d 1340//_____________________________________________________________________________
1341Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1342{
1343 // Reads the trigger decision which is
1344 // stored in Trigger.root file and fills
1345 // the corresponding esd entries
1346
1347 AliInfo("Filling trigger information into the ESD");
1348
1349 if (fRawReader) {
1350 AliCTPRawStream input(fRawReader);
1351 if (!input.Next()) {
1352 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1353 return kFALSE;
1354 }
1355 esd->SetTriggerMask(input.GetClassMask());
1356 esd->SetTriggerCluster(input.GetClusterMask());
1357 }
1358 else {
1359 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1360 if (runloader) {
1361 if (!runloader->LoadTrigger()) {
1362 AliCentralTrigger *aCTP = runloader->GetTrigger();
1363 esd->SetTriggerMask(aCTP->GetClassMask());
1364 esd->SetTriggerCluster(aCTP->GetClusterMask());
1365 }
1366 else {
1367 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1368 return kFALSE;
1369 }
1370 }
1371 else {
1372 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1373 return kFALSE;
1374 }
1375 }
1376
1377 return kTRUE;
1378}
596a855f 1379
001397cd 1380
1381
1382
1383
1384//_____________________________________________________________________________
1385Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1386{
1387 //
1388 // Filling information from RawReader Header
1389 //
1390
1391 AliInfo("Filling information from RawReader Header");
31fd97b2 1392 esd->SetBunchCrossNumber(0);
1393 esd->SetOrbitNumber(0);
001397cd 1394 esd->SetTimeStamp(0);
1395 esd->SetEventType(0);
1396 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1397 if (eventHeader){
31fd97b2 1398 esd->SetBunchCrossNumber((eventHeader->GetP("Id")[0]));
1399 esd->SetOrbitNumber((eventHeader->GetP("Id")[1]));
001397cd 1400 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
31fd97b2 1401 esd->SetEventType((eventHeader->Get("Type")));
001397cd 1402 }
1403
1404 return kTRUE;
1405}
1406
1407
596a855f 1408//_____________________________________________________________________________
1409Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1410{
1411// check whether detName is contained in detectors
1412// if yes, it is removed from detectors
1413
1414 // check if all detectors are selected
1415 if ((detectors.CompareTo("ALL") == 0) ||
1416 detectors.BeginsWith("ALL ") ||
1417 detectors.EndsWith(" ALL") ||
1418 detectors.Contains(" ALL ")) {
1419 detectors = "ALL";
1420 return kTRUE;
1421 }
1422
1423 // search for the given detector
1424 Bool_t result = kFALSE;
1425 if ((detectors.CompareTo(detName) == 0) ||
1426 detectors.BeginsWith(detName+" ") ||
1427 detectors.EndsWith(" "+detName) ||
1428 detectors.Contains(" "+detName+" ")) {
1429 detectors.ReplaceAll(detName, "");
1430 result = kTRUE;
1431 }
1432
1433 // clean up the detectors string
1434 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1435 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1436 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1437
1438 return result;
1439}
e583c30d 1440
f08fc9f5 1441//_____________________________________________________________________________
1442Bool_t AliReconstruction::InitRunLoader()
1443{
1444// get or create the run loader
1445
1446 if (gAlice) delete gAlice;
1447 gAlice = NULL;
1448
b26c3770 1449 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1450 // load all base libraries to get the loader classes
1451 TString libs = gSystem->GetLibraries();
1452 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1453 TString detName = fgkDetectorName[iDet];
1454 if (detName == "HLT") continue;
1455 if (libs.Contains("lib" + detName + "base.so")) continue;
1456 gSystem->Load("lib" + detName + "base.so");
1457 }
f08fc9f5 1458 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1459 if (!fRunLoader) {
1460 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1461 CleanUp();
1462 return kFALSE;
1463 }
b26c3770 1464 fRunLoader->CdGAFile();
1465 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1466 if (fRunLoader->LoadgAlice() == 0) {
1467 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1468 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1469 }
f08fc9f5 1470 }
1471 if (!gAlice && !fRawReader) {
1472 AliError(Form("no gAlice object found in file %s",
1473 fGAliceFileName.Data()));
1474 CleanUp();
1475 return kFALSE;
1476 }
1477
1478 } else { // galice.root does not exist
1479 if (!fRawReader) {
1480 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1481 CleanUp();
1482 return kFALSE;
1483 }
1484 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1485 AliConfig::GetDefaultEventFolderName(),
1486 "recreate");
1487 if (!fRunLoader) {
1488 AliError(Form("could not create run loader in file %s",
1489 fGAliceFileName.Data()));
1490 CleanUp();
1491 return kFALSE;
1492 }
1493 fRunLoader->MakeTree("E");
1494 Int_t iEvent = 0;
1495 while (fRawReader->NextEvent()) {
1496 fRunLoader->SetEventNumber(iEvent);
1497 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1498 iEvent, iEvent);
1499 fRunLoader->MakeTree("H");
1500 fRunLoader->TreeE()->Fill();
1501 iEvent++;
1502 }
1503 fRawReader->RewindEvents();
973388c2 1504 if (fNumberOfEventsPerFile > 0)
1505 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1506 else
1507 fRunLoader->SetNumberOfEventsPerFile(iEvent);
f08fc9f5 1508 fRunLoader->WriteHeader("OVERWRITE");
1509 fRunLoader->CdGAFile();
1510 fRunLoader->Write(0, TObject::kOverwrite);
1511// AliTracker::SetFieldMap(???);
1512 }
1513
1514 return kTRUE;
1515}
1516
c757bafd 1517//_____________________________________________________________________________
b8cd5251 1518AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1519{
f08fc9f5 1520// get the reconstructor object and the loader for a detector
c757bafd 1521
b8cd5251 1522 if (fReconstructor[iDet]) return fReconstructor[iDet];
1523
1524 // load the reconstructor object
1525 TPluginManager* pluginManager = gROOT->GetPluginManager();
1526 TString detName = fgkDetectorName[iDet];
1527 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1528 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1529
1530 if (detName == "HLT") {
1531 if (!gROOT->GetClass("AliLevel3")) {
4aa41877 1532 gSystem->Load("libAliHLTSrc.so");
1533 gSystem->Load("libAliHLTMisc.so");
1534 gSystem->Load("libAliHLTHough.so");
1535 gSystem->Load("libAliHLTComp.so");
b8cd5251 1536 }
1537 }
1538
1539 AliReconstructor* reconstructor = NULL;
1540 // first check if a plugin is defined for the reconstructor
1541 TPluginHandler* pluginHandler =
1542 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1543 // if not, add a plugin for it
1544 if (!pluginHandler) {
b8cd5251 1545 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1546 TString libs = gSystem->GetLibraries();
1547 if (libs.Contains("lib" + detName + "base.so") ||
1548 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1549 pluginManager->AddHandler("AliReconstructor", detName,
1550 recName, detName + "rec", recName + "()");
1551 } else {
1552 pluginManager->AddHandler("AliReconstructor", detName,
1553 recName, detName, recName + "()");
c757bafd 1554 }
b8cd5251 1555 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1556 }
1557 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1558 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1559 }
b8cd5251 1560 if (reconstructor) {
1561 TObject* obj = fOptions.FindObject(detName.Data());
1562 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1563 reconstructor->Init(fRunLoader);
b8cd5251 1564 fReconstructor[iDet] = reconstructor;
1565 }
1566
f08fc9f5 1567 // get or create the loader
1568 if (detName != "HLT") {
1569 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1570 if (!fLoader[iDet]) {
1571 AliConfig::Instance()
1572 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1573 detName, detName);
1574 // first check if a plugin is defined for the loader
bb0901a4 1575 pluginHandler =
f08fc9f5 1576 pluginManager->FindHandler("AliLoader", detName);
1577 // if not, add a plugin for it
1578 if (!pluginHandler) {
1579 TString loaderName = "Ali" + detName + "Loader";
1580 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1581 pluginManager->AddHandler("AliLoader", detName,
1582 loaderName, detName + "base",
1583 loaderName + "(const char*, TFolder*)");
1584 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1585 }
1586 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1587 fLoader[iDet] =
1588 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1589 fRunLoader->GetEventFolder());
1590 }
1591 if (!fLoader[iDet]) { // use default loader
1592 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1593 }
1594 if (!fLoader[iDet]) {
1595 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1596 if (fStopOnError) return NULL;
f08fc9f5 1597 } else {
1598 fRunLoader->AddLoader(fLoader[iDet]);
1599 fRunLoader->CdGAFile();
1600 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1601 fRunLoader->Write(0, TObject::kOverwrite);
1602 }
1603 }
1604 }
1605
b8cd5251 1606 return reconstructor;
c757bafd 1607}
1608
2257f27e 1609//_____________________________________________________________________________
1610Bool_t AliReconstruction::CreateVertexer()
1611{
1612// create the vertexer
1613
b8cd5251 1614 fVertexer = NULL;
1615 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1616 if (itsReconstructor) {
b8cd5251 1617 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1618 }
b8cd5251 1619 if (!fVertexer) {
815c2b38 1620 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1621 if (fStopOnError) return kFALSE;
1622 }
1623
1624 return kTRUE;
1625}
1626
24f7a148 1627//_____________________________________________________________________________
b8cd5251 1628Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1629{
f08fc9f5 1630// create the trackers
24f7a148 1631
b8cd5251 1632 TString detStr = detectors;
1633 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1634 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1635 AliReconstructor* reconstructor = GetReconstructor(iDet);
1636 if (!reconstructor) continue;
1637 TString detName = fgkDetectorName[iDet];
1f46a9ae 1638 if (detName == "HLT") {
1639 fRunHLTTracking = kTRUE;
1640 continue;
1641 }
e66fbafb 1642 if (detName == "MUON") {
1643 fRunMuonTracking = kTRUE;
1644 continue;
1645 }
1646
f08fc9f5 1647
1648 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1649 if (!fTracker[iDet] && (iDet < 7)) {
1650 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1651 if (fStopOnError) return kFALSE;
1652 }
1653 }
1654
24f7a148 1655 return kTRUE;
1656}
1657
e583c30d 1658//_____________________________________________________________________________
b26c3770 1659void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1660{
1661// delete trackers and the run loader and close and delete the file
1662
b8cd5251 1663 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1664 delete fReconstructor[iDet];
1665 fReconstructor[iDet] = NULL;
1666 fLoader[iDet] = NULL;
1667 delete fTracker[iDet];
1668 fTracker[iDet] = NULL;
1669 }
1670 delete fVertexer;
1671 fVertexer = NULL;
9178838a 1672 delete fDiamondProfile;
1673 fDiamondProfile = NULL;
e583c30d 1674
1675 delete fRunLoader;
1676 fRunLoader = NULL;
b649205a 1677 delete fRawReader;
1678 fRawReader = NULL;
e583c30d 1679
1680 if (file) {
1681 file->Close();
1682 delete file;
1683 }
b26c3770 1684
1685 if (fileOld) {
1686 fileOld->Close();
1687 delete fileOld;
1688 gSystem->Unlink("AliESDs.old.root");
1689 }
e583c30d 1690}
24f7a148 1691
1692
1693//_____________________________________________________________________________
1694Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1695{
1696// read the ESD event from a file
1697
1698 if (!esd) return kFALSE;
1699 char fileName[256];
1700 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1701 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1702 if (gSystem->AccessPathName(fileName)) return kFALSE;
1703
f3a97c86 1704 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1705 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1706 TFile* file = TFile::Open(fileName);
1707 if (!file || !file->IsOpen()) {
815c2b38 1708 AliError(Form("opening %s failed", fileName));
24f7a148 1709 delete file;
1710 return kFALSE;
1711 }
1712
1713 gROOT->cd();
1714 delete esd;
1715 esd = (AliESD*) file->Get("ESD");
1716 file->Close();
1717 delete file;
1718 return kTRUE;
1719}
1720
1721//_____________________________________________________________________________
1722void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1723{
1724// write the ESD event to a file
1725
1726 if (!esd) return;
1727 char fileName[256];
1728 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1729 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1730
815c2b38 1731 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1732 TFile* file = TFile::Open(fileName, "recreate");
1733 if (!file || !file->IsOpen()) {
815c2b38 1734 AliError(Form("opening %s failed", fileName));
24f7a148 1735 } else {
1736 esd->Write("ESD");
1737 file->Close();
1738 }
1739 delete file;
1740}
f3a97c86 1741
1742
1743
1744
1745//_____________________________________________________________________________
1746void AliReconstruction::CreateTag(TFile* file)
1747{
2bdb9d38 1748 /////////////
1749 //muon code//
1750 ////////////
56982dd1 1751 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1752 //Variables
56982dd1 1753 Double_t fX,fY,fZ ;
1754 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1755 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1756 Int_t fCharge;
1757 TLorentzVector fEPvector;
1758
1759 Float_t fZVertexCut = 10.0;
1760 Float_t fRhoVertexCut = 2.0;
1761
1762 Float_t fLowPtCut = 1.0;
1763 Float_t fHighPtCut = 3.0;
1764 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1765 ////////////
1766
1767 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1768
089bf903 1769 // Creates the tags for all the events in a given ESD file
4302e20f 1770 Int_t ntrack;
089bf903 1771 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1772 Int_t nPos, nNeg, nNeutr;
1773 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1774 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1775 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1776 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1777 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1778 Int_t fVertexflag;
1387d165 1779 Int_t iRunNumber = 0;
547d75a6 1780 TString fVertexName("default");
4302e20f 1781
f3a97c86 1782 AliRunTag *tag = new AliRunTag();
f3a97c86 1783 AliEventTag *evTag = new AliEventTag();
1784 TTree ttag("T","A Tree with event tags");
2bdb9d38 1785 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1786 btag->SetCompressionLevel(9);
2bdb9d38 1787
f3a97c86 1788 AliInfo(Form("Creating the tags......."));
1789
1790 if (!file || !file->IsOpen()) {
1791 AliError(Form("opening failed"));
1792 delete file;
1793 return ;
2bdb9d38 1794 }
d71aa170 1795 Int_t lastEvent = 0;
f3a97c86 1796 TTree *t = (TTree*) file->Get("esdTree");
1797 TBranch * b = t->GetBranch("ESD");
1798 AliESD *esd = 0;
1799 b->SetAddress(&esd);
1387d165 1800
d71aa170 1801 b->GetEntry(fFirstEvent);
1387d165 1802 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1803
bb0901a4 1804 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1805 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1806 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 1807 ntrack = 0;
1808 nPos = 0;
1809 nNeg = 0;
1810 nNeutr =0;
1811 nK0s = 0;
1812 nNeutrons = 0;
1813 nPi0s = 0;
1814 nGamas = 0;
1815 nProtons = 0;
1816 nKaons = 0;
1817 nPions = 0;
1818 nMuons = 0;
1819 nElectrons = 0;
1820 nCh1GeV = 0;
1821 nCh3GeV = 0;
1822 nCh10GeV = 0;
1823 nMu1GeV = 0;
1824 nMu3GeV = 0;
1825 nMu10GeV = 0;
1826 nEl1GeV = 0;
1827 nEl3GeV = 0;
1828 nEl10GeV = 0;
1829 maxPt = .0;
1830 meanPt = .0;
1831 totalP = .0;
547d75a6 1832 fVertexflag = 0;
d1a50cb5 1833
2bdb9d38 1834 b->GetEntry(iEventNumber);
1387d165 1835 iRunNumber = esd->GetRunNumber();
1836 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1837 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1838 if (!vertexIn) AliError("ESD has not defined vertex.");
1839 if (vertexIn) fVertexName = vertexIn->GetName();
1840 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1841 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1842 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1843 UInt_t status = esdTrack->GetStatus();
f3a97c86 1844
2bdb9d38 1845 //select only tracks with ITS refit
1846 if ((status&AliESDtrack::kITSrefit)==0) continue;
1847 //select only tracks with TPC refit
1848 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1849
2bdb9d38 1850 //select only tracks with the "combined PID"
1851 if ((status&AliESDtrack::kESDpid)==0) continue;
1852 Double_t p[3];
1853 esdTrack->GetPxPyPz(p);
1854 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1855 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1856 totalP += momentum;
1857 meanPt += fPt;
1858 if(fPt > maxPt) maxPt = fPt;
4302e20f 1859
2bdb9d38 1860 if(esdTrack->GetSign() > 0) {
1861 nPos++;
56982dd1 1862 if(fPt > fLowPtCut) nCh1GeV++;
1863 if(fPt > fHighPtCut) nCh3GeV++;
1864 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1865 }
1866 if(esdTrack->GetSign() < 0) {
1867 nNeg++;
56982dd1 1868 if(fPt > fLowPtCut) nCh1GeV++;
1869 if(fPt > fHighPtCut) nCh3GeV++;
1870 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1871 }
1872 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1873
2bdb9d38 1874 //PID
1875 Double_t prob[5];
1876 esdTrack->GetESDpid(prob);
4302e20f 1877
2bdb9d38 1878 Double_t rcc = 0.0;
1879 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1880 if(rcc == 0.0) continue;
1881 //Bayes' formula
1882 Double_t w[5];
1883 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1884
2bdb9d38 1885 //protons
1886 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1887 //kaons
1888 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1889 //pions
1890 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1891 //electrons
1892 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1893 nElectrons++;
56982dd1 1894 if(fPt > fLowPtCut) nEl1GeV++;
1895 if(fPt > fHighPtCut) nEl3GeV++;
1896 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1897 }
1898 ntrack++;
1899 }//track loop
1900
1901 /////////////
1902 //muon code//
1903 ////////////
1904 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1905 // loop over all reconstructed tracks (also first track of combination)
1906 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1907 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1908 if (muonTrack == 0x0) continue;
4302e20f 1909
2bdb9d38 1910 // Coordinates at vertex
56982dd1 1911 fZ = muonTrack->GetZ();
1912 fY = muonTrack->GetBendingCoor();
1913 fX = muonTrack->GetNonBendingCoor();
4302e20f 1914
56982dd1 1915 fThetaX = muonTrack->GetThetaX();
1916 fThetaY = muonTrack->GetThetaY();
4302e20f 1917
56982dd1 1918 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1919 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1920 fPxRec = fPzRec * TMath::Tan(fThetaX);
1921 fPyRec = fPzRec * TMath::Tan(fThetaY);
1922 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1923
2bdb9d38 1924 //ChiSquare of the track if needed
56982dd1 1925 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1926 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1927 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1928
2bdb9d38 1929 // total number of muons inside a vertex cut
56982dd1 1930 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1931 nMuons++;
56982dd1 1932 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1933 nMu1GeV++;
56982dd1 1934 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1935 nMu3GeV++;
56982dd1 1936 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1937 nMu10GeV++;
1938 }
1939 }
1940 }
1941 }
1942 }//muon track loop
1943
1944 // Fill the event tags
1945 if(ntrack != 0)
1946 meanPt = meanPt/ntrack;
1947
1948 evTag->SetEventId(iEventNumber+1);
547d75a6 1949 if (vertexIn) {
1950 evTag->SetVertexX(vertexIn->GetXv());
1951 evTag->SetVertexY(vertexIn->GetYv());
1952 evTag->SetVertexZ(vertexIn->GetZv());
1953 evTag->SetVertexZError(vertexIn->GetZRes());
1954 }
d1a50cb5 1955 evTag->SetVertexFlag(fVertexflag);
1956
2bdb9d38 1957 evTag->SetT0VertexZ(esd->GetT0zVertex());
1958
8bd8ac26 1959 evTag->SetTriggerMask(esd->GetTriggerMask());
1960 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1961
32a5cab4 1962 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1963 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1964 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1965 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1966 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1967 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1968
1969
1970 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1971 evTag->SetNumOfPosTracks(nPos);
1972 evTag->SetNumOfNegTracks(nNeg);
1973 evTag->SetNumOfNeutrTracks(nNeutr);
1974
1975 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1976 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1977 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1978 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1979
1980 evTag->SetNumOfProtons(nProtons);
1981 evTag->SetNumOfKaons(nKaons);
1982 evTag->SetNumOfPions(nPions);
1983 evTag->SetNumOfMuons(nMuons);
1984 evTag->SetNumOfElectrons(nElectrons);
1985 evTag->SetNumOfPhotons(nGamas);
1986 evTag->SetNumOfPi0s(nPi0s);
1987 evTag->SetNumOfNeutrons(nNeutrons);
1988 evTag->SetNumOfKaon0s(nK0s);
1989
1990 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1991 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1992 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1993 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1994 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1995 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1996 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1997 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1998 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1999
85c60a8e 2000 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
2001 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 2002
2003 evTag->SetTotalMomentum(totalP);
2004 evTag->SetMeanPt(meanPt);
2005 evTag->SetMaxPt(maxPt);
2006
1387d165 2007 tag->SetRunId(iInitRunNumber);
2bdb9d38 2008 tag->AddEventTag(*evTag);
2009 }
bb0901a4 2010 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 2011 else lastEvent = fLastEvent;
f3a97c86 2012
2013 ttag.Fill();
2014 tag->Clear();
2015
2016 char fileName[256];
2017 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 2018 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 2019 AliInfo(Form("writing tags to file %s", fileName));
2020 AliDebug(1, Form("writing tags to file %s", fileName));
2021
2022 TFile* ftag = TFile::Open(fileName, "recreate");
2023 ftag->cd();
2024 ttag.Write();
2025 ftag->Close();
2026 file->cd();
2027 delete tag;
f3a97c86 2028 delete evTag;
2029}
2030
a7807689 2031//_____________________________________________________________________________
2032void AliReconstruction::CreateAOD(TFile* esdFile)
2033{
2034 // do nothing for now
2035
2036 return;
2037}
2038
2039
98937d93 2040void AliReconstruction::WriteAlignmentData(AliESD* esd)
2041{
2042 // Write space-points which are then used in the alignment procedures
2043 // For the moment only ITS, TRD and TPC
2044
2045 // Load TOF clusters
d528ee75 2046 if (fTracker[3]){
2047 fLoader[3]->LoadRecPoints("read");
2048 TTree* tree = fLoader[3]->TreeR();
2049 if (!tree) {
2050 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2051 return;
2052 }
2053 fTracker[3]->LoadClusters(tree);
98937d93 2054 }
98937d93 2055 Int_t ntracks = esd->GetNumberOfTracks();
2056 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2057 {
2058 AliESDtrack *track = esd->GetTrack(itrack);
2059 Int_t nsp = 0;
ef7253ac 2060 Int_t idx[200];
98937d93 2061 for (Int_t iDet = 3; iDet >= 0; iDet--)
2062 nsp += track->GetNcls(iDet);
2063 if (nsp) {
2064 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2065 track->SetTrackPointArray(sp);
2066 Int_t isptrack = 0;
2067 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2068 AliTracker *tracker = fTracker[iDet];
2069 if (!tracker) continue;
2070 Int_t nspdet = track->GetNcls(iDet);
98937d93 2071 if (nspdet <= 0) continue;
2072 track->GetClusters(iDet,idx);
2073 AliTrackPoint p;
2074 Int_t isp = 0;
2075 Int_t isp2 = 0;
2076 while (isp < nspdet) {
2077 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 2078 const Int_t kNTPCmax = 159;
2079 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2080 if (!isvalid) continue;
2081 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2082 }
98937d93 2083 }
2084 }
2085 }
d528ee75 2086 if (fTracker[3]){
2087 fTracker[3]->UnloadClusters();
2088 fLoader[3]->UnloadRecPoints();
2089 }
98937d93 2090}
2e3550da 2091
2092//_____________________________________________________________________________
2093void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2094{
2095 // The method reads the raw-data error log
2096 // accumulated within the rawReader.
2097 // It extracts the raw-data errors related to
2098 // the current event and stores them into
2099 // a TClonesArray inside the esd object.
2100
2101 if (!fRawReader) return;
2102
2103 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2104
2105 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2106 if (!log) continue;
2107 if (iEvent != log->GetEventNumber()) continue;
2108
2109 esd->AddRawDataErrorLog(log);
2110 }
2111
2112}