]> 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
A macro to create a nominal beam position object in OCDB
[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;
a03b0371 539
540 // Update the TGeoPhysicalNodes
541 gGeoManager->RefreshPhysicalNodes();
542
6bae477a 543 return kTRUE;
544}
596a855f 545
546//_____________________________________________________________________________
547void AliReconstruction::SetGAliceFile(const char* fileName)
548{
549// set the name of the galice file
550
551 fGAliceFileName = fileName;
552}
553
efd2085e 554//_____________________________________________________________________________
555void AliReconstruction::SetOption(const char* detector, const char* option)
556{
557// set options for the reconstruction of a detector
558
559 TObject* obj = fOptions.FindObject(detector);
560 if (obj) fOptions.Remove(obj);
561 fOptions.Add(new TNamed(detector, option));
562}
563
596a855f 564
565//_____________________________________________________________________________
4a33489c 566Bool_t AliReconstruction::Run(const char* input)
596a855f 567{
568// run the reconstruction
569
b649205a 570 // set the input
571 if (!input) input = fInput.Data();
572 TString fileName(input);
573 if (fileName.EndsWith("/")) {
574 fRawReader = new AliRawReaderFile(fileName);
575 } else if (fileName.EndsWith(".root")) {
576 fRawReader = new AliRawReaderRoot(fileName);
577 } else if (!fileName.IsNull()) {
578 fRawReader = new AliRawReaderDate(fileName);
579 fRawReader->SelectEvents(7);
580 }
35042093 581 if (!fEquipIdMap.IsNull() && fRawReader)
582 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
583
b649205a 584
f08fc9f5 585 // get the run loader
586 if (!InitRunLoader()) return kFALSE;
596a855f 587
ec92bee0 588 // Initialize the CDB storage
589 InitCDBStorage();
590
6bae477a 591 // Set run number in CDBManager (if it is not already set by the user)
592 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
593
594 // Import ideal TGeo geometry and apply misalignment
595 if (!gGeoManager) {
596 TString geom(gSystem->DirName(fGAliceFileName));
597 geom += "/geometry.root";
598 TGeoManager::Import(geom.Data());
599 if (!gGeoManager) if (fStopOnError) return kFALSE;
600 }
8e245d15 601
602 AliCDBManager* man = AliCDBManager::Instance();
6bae477a 603 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
604
596a855f 605 // local reconstruction
59697224 606 if (!fRunLocalReconstruction.IsNull()) {
607 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 608 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 609 }
610 }
b26c3770 611// if (!fRunVertexFinder && fRunTracking.IsNull() &&
612// fFillESD.IsNull()) return kTRUE;
2257f27e 613
614 // get vertexer
615 if (fRunVertexFinder && !CreateVertexer()) {
616 if (fStopOnError) {
617 CleanUp();
618 return kFALSE;
619 }
620 }
596a855f 621
f08fc9f5 622 // get trackers
b8cd5251 623 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 624 if (fStopOnError) {
625 CleanUp();
626 return kFALSE;
627 }
596a855f 628 }
24f7a148 629
9db3a215 630
631 TStopwatch stopwatch;
632 stopwatch.Start();
633
b26c3770 634 // get the possibly already existing ESD file and tree
1f46a9ae 635 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
b26c3770 636 TFile* fileOld = NULL;
1f46a9ae 637 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 638 if (!gSystem->AccessPathName("AliESDs.root")){
639 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
640 fileOld = TFile::Open("AliESDs.old.root");
641 if (fileOld && fileOld->IsOpen()) {
642 treeOld = (TTree*) fileOld->Get("esdTree");
643 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
1f46a9ae 644 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
645 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
b26c3770 646 }
647 }
648
36711aa4 649 // create the ESD output file and tree
596a855f 650 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
651 if (!file->IsOpen()) {
815c2b38 652 AliError("opening AliESDs.root failed");
b26c3770 653 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 654 }
36711aa4 655 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
656 tree->Branch("ESD", "AliESD", &esd);
1f46a9ae 657 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
658 hlttree->Branch("ESD", "AliESD", &hltesd);
659 delete esd; delete hltesd;
660 esd = NULL; hltesd = NULL;
1d99986f 661
500d54ab 662 // create the branch with ESD additions
663 AliESDfriend *esdf=0;
1d99986f 664 if (fWriteESDfriend) {
500d54ab 665 TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
666 br->SetFile("AliESDfriends.root");
1d99986f 667 }
668
20e5681c 669 AliVertexerTracks tVertexer(AliTracker::GetBz());
9178838a 670 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
c5e3e5d1 671
596a855f 672 // loop over events
b649205a 673 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 674
596a855f 675 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 676 if (fRawReader) fRawReader->NextEvent();
4a33489c 677 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
b26c3770 678 // copy old ESD to the new one
679 if (treeOld) {
680 treeOld->SetBranchAddress("ESD", &esd);
681 treeOld->GetEntry(iEvent);
682 }
683 tree->Fill();
1f46a9ae 684 if (hlttreeOld) {
685 hlttreeOld->SetBranchAddress("ESD", &hltesd);
686 hlttreeOld->GetEntry(iEvent);
687 }
688 hlttree->Fill();
b26c3770 689 continue;
690 }
691
815c2b38 692 AliInfo(Form("processing event %d", iEvent));
596a855f 693 fRunLoader->GetEvent(iEvent);
24f7a148 694
bb0901a4 695 char aFileName[256];
696 sprintf(aFileName, "ESD_%d.%d_final.root",
f08fc9f5 697 fRunLoader->GetHeader()->GetRun(),
698 fRunLoader->GetHeader()->GetEventNrInRun());
bb0901a4 699 if (!gSystem->AccessPathName(aFileName)) continue;
24f7a148 700
b26c3770 701 // local reconstruction
702 if (!fRunLocalReconstruction.IsNull()) {
703 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
704 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
705 }
706 }
707
1f46a9ae 708 esd = new AliESD; hltesd = new AliESD;
f08fc9f5 709 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 710 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
31fd97b2 711 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
712 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
d6ee376f 713
714 // Set magnetic field from the tracker
715 esd->SetMagneticField(AliTracker::GetBz());
716 hltesd->SetMagneticField(AliTracker::GetBz());
596a855f 717
2e3550da 718 // Fill raw-data error log into the ESD
719 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
720
2257f27e 721 // vertex finder
722 if (fRunVertexFinder) {
723 if (!ReadESD(esd, "vertex")) {
724 if (!RunVertexFinder(esd)) {
b26c3770 725 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 726 }
727 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
728 }
729 }
730
1f46a9ae 731 // HLT tracking
732 if (!fRunTracking.IsNull()) {
733 if (fRunHLTTracking) {
734 hltesd->SetVertex(esd->GetVertex());
735 if (!RunHLTTracking(hltesd)) {
736 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
737 }
738 }
739 }
740
e66fbafb 741 // Muon tracking
b8cd5251 742 if (!fRunTracking.IsNull()) {
e66fbafb 743 if (fRunMuonTracking) {
744 if (!RunMuonTracking()) {
b26c3770 745 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 746 }
596a855f 747 }
748 }
749
e66fbafb 750 // barrel tracking
751 if (!fRunTracking.IsNull()) {
21c573b7 752 if (!ReadESD(esd, "tracking")) {
753 if (!RunTracking(esd)) {
754 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
e66fbafb 755 }
21c573b7 756 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
e66fbafb 757 }
758 }
21c573b7 759
596a855f 760 // fill ESD
761 if (!fFillESD.IsNull()) {
762 if (!FillESD(esd, fFillESD)) {
b26c3770 763 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 764 }
765 }
001397cd 766 // fill Event header information from the RawEventHeader
767 if (fRawReader){FillRawEventHeaderESD(esd);}
596a855f 768
769 // combined PID
770 AliESDpid::MakePID(esd);
24f7a148 771 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 772
b647652d 773 if (fFillTriggerESD) {
774 if (!ReadESD(esd, "trigger")) {
775 if (!FillTriggerESD(esd)) {
776 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
777 }
778 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
779 }
780 }
781
a03cd2e0 782 //Try to improve the reconstructed primary vertex position using the tracks
783 AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
784 if (pvtx)
785 if (pvtx->GetStatus()) {
786 // Store the improved primary vertex
787 esd->SetPrimaryVertex(pvtx);
788 // Propagate the tracks to the DCA to the improved primary vertex
789 Double_t somethingbig = 777.;
790 Double_t bz = esd->GetMagneticField();
791 Int_t nt=esd->GetNumberOfTracks();
792 while (nt--) {
793 AliESDtrack *t = esd->GetTrack(nt);
794 t->RelateToVertex(pvtx, bz, somethingbig);
795 }
796 }
c5e3e5d1 797
5e4ff34d 798 {
799 // V0 finding
800 AliV0vertexer vtxer;
801 vtxer.Tracks2V0vertices(esd);
802
803 // Cascade finding
804 AliCascadeVertexer cvtxer;
805 cvtxer.V0sTracks2CascadeVertices(esd);
806 }
807
596a855f 808 // write ESD
1d99986f 809 if (fWriteESDfriend) {
d75007f6 810 esdf=new AliESDfriend();
811 esd->GetESDfriend(esdf);
1d99986f 812 }
500d54ab 813 tree->Fill();
814
815 // write HLT ESD
816 hlttree->Fill();
1d99986f 817
f3a97c86 818 if (fCheckPointLevel > 0) WriteESD(esd, "final");
819
500d54ab 820 delete esd; delete esdf; delete hltesd;
821 esd = NULL; esdf=NULL; hltesd = NULL;
596a855f 822 }
823
9db3a215 824 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
825 stopwatch.RealTime(),stopwatch.CpuTime()));
826
36711aa4 827 file->cd();
a9c0e6db 828 if (fWriteESDfriend)
829 tree->SetBranchStatus("ESDfriend*",0);
36711aa4 830 tree->Write();
1f46a9ae 831 hlttree->Write();
f3a97c86 832
a7807689 833 if (fWriteAOD) {
834 CreateAOD(file);
835 }
836
f3a97c86 837 // Create tags for the events in the ESD tree (the ESD tree is always present)
838 // In case of empty events the tags will contain dummy values
839 CreateTag(file);
b26c3770 840 CleanUp(file, fileOld);
596a855f 841
842 return kTRUE;
843}
844
845
846//_____________________________________________________________________________
59697224 847Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 848{
59697224 849// run the local reconstruction
596a855f 850
030b532d 851 TStopwatch stopwatch;
852 stopwatch.Start();
853
8e245d15 854 AliCDBManager* man = AliCDBManager::Instance();
855 Bool_t origCache = man->GetCacheFlag();
856
596a855f 857 TString detStr = detectors;
b8cd5251 858 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
859 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
860 AliReconstructor* reconstructor = GetReconstructor(iDet);
861 if (!reconstructor) continue;
b26c3770 862 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 863
864 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
865 TStopwatch stopwatchDet;
866 stopwatchDet.Start();
8e245d15 867
868 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
869
870 man->SetCacheFlag(kTRUE);
871 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
872 man->GetAll(calibPath); // entries are cached!
873
b8cd5251 874 if (fRawReader) {
875 fRawReader->RewindEvents();
876 reconstructor->Reconstruct(fRunLoader, fRawReader);
877 } else {
878 reconstructor->Reconstruct(fRunLoader);
596a855f 879 }
5f8272e1 880 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
881 fgkDetectorName[iDet],
882 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
8e245d15 883
884 // unload calibration data
885 man->ClearCache();
596a855f 886 }
887
8e245d15 888 man->SetCacheFlag(origCache);
889
596a855f 890 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 891 AliError(Form("the following detectors were not found: %s",
892 detStr.Data()));
596a855f 893 if (fStopOnError) return kFALSE;
894 }
895
5f8272e1 896 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
897 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 898
596a855f 899 return kTRUE;
900}
901
b26c3770 902//_____________________________________________________________________________
903Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
904{
905// run the local reconstruction
906
907 TStopwatch stopwatch;
908 stopwatch.Start();
909
910 TString detStr = detectors;
911 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
912 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
913 AliReconstructor* reconstructor = GetReconstructor(iDet);
914 if (!reconstructor) continue;
915 AliLoader* loader = fLoader[iDet];
916
917 // conversion of digits
918 if (fRawReader && reconstructor->HasDigitConversion()) {
919 AliInfo(Form("converting raw data digits into root objects for %s",
920 fgkDetectorName[iDet]));
921 TStopwatch stopwatchDet;
922 stopwatchDet.Start();
923 loader->LoadDigits("update");
924 loader->CleanDigits();
925 loader->MakeDigitsContainer();
926 TTree* digitsTree = loader->TreeD();
927 reconstructor->ConvertDigits(fRawReader, digitsTree);
928 loader->WriteDigits("OVERWRITE");
929 loader->UnloadDigits();
5f8272e1 930 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
931 fgkDetectorName[iDet],
932 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 933 }
934
935 // local reconstruction
936 if (!reconstructor->HasLocalReconstruction()) continue;
937 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
938 TStopwatch stopwatchDet;
939 stopwatchDet.Start();
940 loader->LoadRecPoints("update");
941 loader->CleanRecPoints();
942 loader->MakeRecPointsContainer();
943 TTree* clustersTree = loader->TreeR();
944 if (fRawReader && !reconstructor->HasDigitConversion()) {
945 reconstructor->Reconstruct(fRawReader, clustersTree);
946 } else {
947 loader->LoadDigits("read");
948 TTree* digitsTree = loader->TreeD();
949 if (!digitsTree) {
950 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
951 if (fStopOnError) return kFALSE;
952 } else {
953 reconstructor->Reconstruct(digitsTree, clustersTree);
954 }
955 loader->UnloadDigits();
956 }
957 loader->WriteRecPoints("OVERWRITE");
958 loader->UnloadRecPoints();
5f8272e1 959 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
960 fgkDetectorName[iDet],
961 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 962 }
963
964 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
965 AliError(Form("the following detectors were not found: %s",
966 detStr.Data()));
967 if (fStopOnError) return kFALSE;
968 }
5f8272e1 969
970 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
971 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 972
973 return kTRUE;
974}
975
596a855f 976//_____________________________________________________________________________
2257f27e 977Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 978{
979// run the barrel tracking
980
030b532d 981 TStopwatch stopwatch;
982 stopwatch.Start();
983
2257f27e 984 AliESDVertex* vertex = NULL;
985 Double_t vtxPos[3] = {0, 0, 0};
986 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
987 TArrayF mcVertex(3);
a6b0b91b 988 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
989 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
990 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
991 }
2257f27e 992
b8cd5251 993 if (fVertexer) {
815c2b38 994 AliInfo("running the ITS vertex finder");
b26c3770 995 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 996 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 997 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 998 if(!vertex){
815c2b38 999 AliWarning("Vertex not found");
c710f220 1000 vertex = new AliESDVertex();
d1a50cb5 1001 vertex->SetName("default");
2257f27e 1002 }
1003 else {
1004 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 1005 vertex->SetName("reconstructed");
2257f27e 1006 }
1007
1008 } else {
815c2b38 1009 AliInfo("getting the primary vertex from MC");
2257f27e 1010 vertex = new AliESDVertex(vtxPos, vtxErr);
1011 }
1012
1013 if (vertex) {
1014 vertex->GetXYZ(vtxPos);
1015 vertex->GetSigmaXYZ(vtxErr);
1016 } else {
815c2b38 1017 AliWarning("no vertex reconstructed");
2257f27e 1018 vertex = new AliESDVertex(vtxPos, vtxErr);
1019 }
1020 esd->SetVertex(vertex);
32e449be 1021 // if SPD multiplicity has been determined, it is stored in the ESD
1022 AliMultiplicity *mult= fVertexer->GetMultiplicity();
1023 if(mult)esd->SetMultiplicity(mult);
1024
b8cd5251 1025 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1026 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1027 }
2257f27e 1028 delete vertex;
1029
5f8272e1 1030 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1031 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 1032
1033 return kTRUE;
1034}
1035
1f46a9ae 1036//_____________________________________________________________________________
1037Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1038{
1039// run the HLT barrel tracking
1040
1041 TStopwatch stopwatch;
1042 stopwatch.Start();
1043
1044 if (!fRunLoader) {
1045 AliError("Missing runLoader!");
1046 return kFALSE;
1047 }
1048
1049 AliInfo("running HLT tracking");
1050
1051 // Get a pointer to the HLT reconstructor
1052 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1053 if (!reconstructor) return kFALSE;
1054
1055 // TPC + ITS
1056 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1057 TString detName = fgkDetectorName[iDet];
1058 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1059 reconstructor->SetOption(detName.Data());
1060 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1061 if (!tracker) {
1062 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1063 if (fStopOnError) return kFALSE;
9dcc06e1 1064 continue;
1f46a9ae 1065 }
1066 Double_t vtxPos[3];
1067 Double_t vtxErr[3]={0.005,0.005,0.010};
1068 const AliESDVertex *vertex = esd->GetVertex();
1069 vertex->GetXYZ(vtxPos);
1070 tracker->SetVertex(vtxPos,vtxErr);
1071 if(iDet != 1) {
1072 fLoader[iDet]->LoadRecPoints("read");
1073 TTree* tree = fLoader[iDet]->TreeR();
1074 if (!tree) {
1075 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1076 return kFALSE;
1077 }
1078 tracker->LoadClusters(tree);
1079 }
1080 if (tracker->Clusters2Tracks(esd) != 0) {
1081 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1082 return kFALSE;
1083 }
1084 if(iDet != 1) {
1085 tracker->UnloadClusters();
1086 }
1087 delete tracker;
1088 }
1089
5f8272e1 1090 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1091 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 1092
1093 return kTRUE;
1094}
1095
e66fbafb 1096//_____________________________________________________________________________
1097Bool_t AliReconstruction::RunMuonTracking()
1098{
1099// run the muon spectrometer tracking
1100
1101 TStopwatch stopwatch;
1102 stopwatch.Start();
1103
1104 if (!fRunLoader) {
1105 AliError("Missing runLoader!");
1106 return kFALSE;
1107 }
1108 Int_t iDet = 7; // for MUON
1109
1110 AliInfo("is running...");
1111
1112 // Get a pointer to the MUON reconstructor
1113 AliReconstructor *reconstructor = GetReconstructor(iDet);
1114 if (!reconstructor) return kFALSE;
1115
1116
1117 TString detName = fgkDetectorName[iDet];
1118 AliDebug(1, Form("%s tracking", detName.Data()));
1119 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1120 if (!tracker) {
1121 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1122 return kFALSE;
1123 }
1124
1125 // create Tracks
1126 fLoader[iDet]->LoadTracks("update");
1127 fLoader[iDet]->CleanTracks();
1128 fLoader[iDet]->MakeTracksContainer();
1129
1130 // read RecPoints
1131 fLoader[iDet]->LoadRecPoints("read");
1132
1133 if (!tracker->Clusters2Tracks(0x0)) {
1134 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1135 return kFALSE;
1136 }
1137 fLoader[iDet]->UnloadRecPoints();
1138
1139 fLoader[iDet]->WriteTracks("OVERWRITE");
1140 fLoader[iDet]->UnloadTracks();
1141
1142 delete tracker;
1143
1144
1145 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1146 stopwatch.RealTime(),stopwatch.CpuTime()));
1147
1148 return kTRUE;
1149}
1150
1151
2257f27e 1152//_____________________________________________________________________________
1153Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1154{
1155// run the barrel tracking
1156
1157 TStopwatch stopwatch;
1158 stopwatch.Start();
24f7a148 1159
815c2b38 1160 AliInfo("running tracking");
596a855f 1161
b8cd5251 1162 // pass 1: TPC + ITS inwards
1163 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1164 if (!fTracker[iDet]) continue;
1165 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 1166
b8cd5251 1167 // load clusters
1168 fLoader[iDet]->LoadRecPoints("read");
1169 TTree* tree = fLoader[iDet]->TreeR();
1170 if (!tree) {
1171 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1172 return kFALSE;
1173 }
b8cd5251 1174 fTracker[iDet]->LoadClusters(tree);
1175
1176 // run tracking
1177 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1178 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 1179 return kFALSE;
1180 }
b8cd5251 1181 if (fCheckPointLevel > 1) {
1182 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1183 }
878e1fe1 1184 // preliminary PID in TPC needed by the ITS tracker
1185 if (iDet == 1) {
1186 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 1187 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 1188 AliESDpid::MakePID(esd);
1189 }
b8cd5251 1190 }
596a855f 1191
b8cd5251 1192 // pass 2: ALL backwards
1193 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1194 if (!fTracker[iDet]) continue;
1195 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1196
1197 // load clusters
1198 if (iDet > 1) { // all except ITS, TPC
1199 TTree* tree = NULL;
7b61cd9c 1200 fLoader[iDet]->LoadRecPoints("read");
1201 tree = fLoader[iDet]->TreeR();
b8cd5251 1202 if (!tree) {
1203 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1204 return kFALSE;
1205 }
b8cd5251 1206 fTracker[iDet]->LoadClusters(tree);
1207 }
24f7a148 1208
b8cd5251 1209 // run tracking
1210 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1211 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1212 return kFALSE;
1213 }
1214 if (fCheckPointLevel > 1) {
1215 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1216 }
24f7a148 1217
b8cd5251 1218 // unload clusters
1219 if (iDet > 2) { // all except ITS, TPC, TRD
1220 fTracker[iDet]->UnloadClusters();
7b61cd9c 1221 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1222 }
8f37df88 1223 // updated PID in TPC needed by the ITS tracker -MI
1224 if (iDet == 1) {
1225 GetReconstructor(1)->FillESD(fRunLoader, esd);
1226 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1227 AliESDpid::MakePID(esd);
1228 }
b8cd5251 1229 }
596a855f 1230
98937d93 1231 // write space-points to the ESD in case alignment data output
1232 // is switched on
1233 if (fWriteAlignmentData)
1234 WriteAlignmentData(esd);
1235
b8cd5251 1236 // pass 3: TRD + TPC + ITS refit inwards
1237 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1238 if (!fTracker[iDet]) continue;
1239 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 1240
b8cd5251 1241 // run tracking
1242 if (fTracker[iDet]->RefitInward(esd) != 0) {
1243 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1244 return kFALSE;
1245 }
1246 if (fCheckPointLevel > 1) {
1247 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1248 }
596a855f 1249
b8cd5251 1250 // unload clusters
1251 fTracker[iDet]->UnloadClusters();
1252 fLoader[iDet]->UnloadRecPoints();
1253 }
ff8bb5ae 1254 //
1255 // Propagate track to the vertex - if not done by ITS
1256 //
1257 Int_t ntracks = esd->GetNumberOfTracks();
1258 for (Int_t itrack=0; itrack<ntracks; itrack++){
1259 const Double_t kRadius = 3; // beam pipe radius
1260 const Double_t kMaxStep = 5; // max step
1261 const Double_t kMaxD = 123456; // max distance to prim vertex
1262 Double_t fieldZ = AliTracker::GetBz(); //
1263 AliESDtrack * track = esd->GetTrack(itrack);
1264 if (!track) continue;
1265 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
a28195f7 1266 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
ff8bb5ae 1267 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1268 }
1269
5f8272e1 1270 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1271 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1272
596a855f 1273 return kTRUE;
1274}
1275
1276//_____________________________________________________________________________
24f7a148 1277Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1278{
1279// fill the event summary data
1280
030b532d 1281 TStopwatch stopwatch;
1282 stopwatch.Start();
815c2b38 1283 AliInfo("filling ESD");
030b532d 1284
596a855f 1285 TString detStr = detectors;
b8cd5251 1286 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1287 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1288 AliReconstructor* reconstructor = GetReconstructor(iDet);
1289 if (!reconstructor) continue;
1290
1291 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1292 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 1293 TTree* clustersTree = NULL;
1294 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1295 fLoader[iDet]->LoadRecPoints("read");
1296 clustersTree = fLoader[iDet]->TreeR();
1297 if (!clustersTree) {
1298 AliError(Form("Can't get the %s clusters tree",
1299 fgkDetectorName[iDet]));
1300 if (fStopOnError) return kFALSE;
1301 }
1302 }
1303 if (fRawReader && !reconstructor->HasDigitConversion()) {
1304 reconstructor->FillESD(fRawReader, clustersTree, esd);
1305 } else {
1306 TTree* digitsTree = NULL;
1307 if (fLoader[iDet]) {
1308 fLoader[iDet]->LoadDigits("read");
1309 digitsTree = fLoader[iDet]->TreeD();
1310 if (!digitsTree) {
1311 AliError(Form("Can't get the %s digits tree",
1312 fgkDetectorName[iDet]));
1313 if (fStopOnError) return kFALSE;
1314 }
1315 }
1316 reconstructor->FillESD(digitsTree, clustersTree, esd);
1317 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1318 }
1319 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1320 fLoader[iDet]->UnloadRecPoints();
1321 }
1322
b8cd5251 1323 if (fRawReader) {
1324 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1325 } else {
1326 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1327 }
b8cd5251 1328 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1329 }
1330 }
1331
1332 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1333 AliError(Form("the following detectors were not found: %s",
1334 detStr.Data()));
596a855f 1335 if (fStopOnError) return kFALSE;
1336 }
1337
5f8272e1 1338 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1339 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1340
596a855f 1341 return kTRUE;
1342}
1343
b647652d 1344//_____________________________________________________________________________
1345Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1346{
1347 // Reads the trigger decision which is
1348 // stored in Trigger.root file and fills
1349 // the corresponding esd entries
1350
1351 AliInfo("Filling trigger information into the ESD");
1352
1353 if (fRawReader) {
1354 AliCTPRawStream input(fRawReader);
1355 if (!input.Next()) {
1356 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1357 return kFALSE;
1358 }
1359 esd->SetTriggerMask(input.GetClassMask());
1360 esd->SetTriggerCluster(input.GetClusterMask());
1361 }
1362 else {
1363 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1364 if (runloader) {
1365 if (!runloader->LoadTrigger()) {
1366 AliCentralTrigger *aCTP = runloader->GetTrigger();
1367 esd->SetTriggerMask(aCTP->GetClassMask());
1368 esd->SetTriggerCluster(aCTP->GetClusterMask());
1369 }
1370 else {
1371 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1372 return kFALSE;
1373 }
1374 }
1375 else {
1376 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1377 return kFALSE;
1378 }
1379 }
1380
1381 return kTRUE;
1382}
596a855f 1383
001397cd 1384
1385
1386
1387
1388//_____________________________________________________________________________
1389Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1390{
1391 //
1392 // Filling information from RawReader Header
1393 //
1394
1395 AliInfo("Filling information from RawReader Header");
31fd97b2 1396 esd->SetBunchCrossNumber(0);
1397 esd->SetOrbitNumber(0);
001397cd 1398 esd->SetTimeStamp(0);
1399 esd->SetEventType(0);
1400 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1401 if (eventHeader){
31fd97b2 1402 esd->SetBunchCrossNumber((eventHeader->GetP("Id")[0]));
1403 esd->SetOrbitNumber((eventHeader->GetP("Id")[1]));
001397cd 1404 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
31fd97b2 1405 esd->SetEventType((eventHeader->Get("Type")));
001397cd 1406 }
1407
1408 return kTRUE;
1409}
1410
1411
596a855f 1412//_____________________________________________________________________________
1413Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1414{
1415// check whether detName is contained in detectors
1416// if yes, it is removed from detectors
1417
1418 // check if all detectors are selected
1419 if ((detectors.CompareTo("ALL") == 0) ||
1420 detectors.BeginsWith("ALL ") ||
1421 detectors.EndsWith(" ALL") ||
1422 detectors.Contains(" ALL ")) {
1423 detectors = "ALL";
1424 return kTRUE;
1425 }
1426
1427 // search for the given detector
1428 Bool_t result = kFALSE;
1429 if ((detectors.CompareTo(detName) == 0) ||
1430 detectors.BeginsWith(detName+" ") ||
1431 detectors.EndsWith(" "+detName) ||
1432 detectors.Contains(" "+detName+" ")) {
1433 detectors.ReplaceAll(detName, "");
1434 result = kTRUE;
1435 }
1436
1437 // clean up the detectors string
1438 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1439 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1440 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1441
1442 return result;
1443}
e583c30d 1444
f08fc9f5 1445//_____________________________________________________________________________
1446Bool_t AliReconstruction::InitRunLoader()
1447{
1448// get or create the run loader
1449
1450 if (gAlice) delete gAlice;
1451 gAlice = NULL;
1452
b26c3770 1453 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1454 // load all base libraries to get the loader classes
1455 TString libs = gSystem->GetLibraries();
1456 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1457 TString detName = fgkDetectorName[iDet];
1458 if (detName == "HLT") continue;
1459 if (libs.Contains("lib" + detName + "base.so")) continue;
1460 gSystem->Load("lib" + detName + "base.so");
1461 }
f08fc9f5 1462 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1463 if (!fRunLoader) {
1464 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1465 CleanUp();
1466 return kFALSE;
1467 }
b26c3770 1468 fRunLoader->CdGAFile();
1469 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1470 if (fRunLoader->LoadgAlice() == 0) {
1471 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1472 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1473 }
f08fc9f5 1474 }
1475 if (!gAlice && !fRawReader) {
1476 AliError(Form("no gAlice object found in file %s",
1477 fGAliceFileName.Data()));
1478 CleanUp();
1479 return kFALSE;
1480 }
1481
1482 } else { // galice.root does not exist
1483 if (!fRawReader) {
1484 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1485 CleanUp();
1486 return kFALSE;
1487 }
1488 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1489 AliConfig::GetDefaultEventFolderName(),
1490 "recreate");
1491 if (!fRunLoader) {
1492 AliError(Form("could not create run loader in file %s",
1493 fGAliceFileName.Data()));
1494 CleanUp();
1495 return kFALSE;
1496 }
1497 fRunLoader->MakeTree("E");
1498 Int_t iEvent = 0;
1499 while (fRawReader->NextEvent()) {
1500 fRunLoader->SetEventNumber(iEvent);
1501 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1502 iEvent, iEvent);
1503 fRunLoader->MakeTree("H");
1504 fRunLoader->TreeE()->Fill();
1505 iEvent++;
1506 }
1507 fRawReader->RewindEvents();
973388c2 1508 if (fNumberOfEventsPerFile > 0)
1509 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1510 else
1511 fRunLoader->SetNumberOfEventsPerFile(iEvent);
f08fc9f5 1512 fRunLoader->WriteHeader("OVERWRITE");
1513 fRunLoader->CdGAFile();
1514 fRunLoader->Write(0, TObject::kOverwrite);
1515// AliTracker::SetFieldMap(???);
1516 }
1517
1518 return kTRUE;
1519}
1520
c757bafd 1521//_____________________________________________________________________________
b8cd5251 1522AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1523{
f08fc9f5 1524// get the reconstructor object and the loader for a detector
c757bafd 1525
b8cd5251 1526 if (fReconstructor[iDet]) return fReconstructor[iDet];
1527
1528 // load the reconstructor object
1529 TPluginManager* pluginManager = gROOT->GetPluginManager();
1530 TString detName = fgkDetectorName[iDet];
1531 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1532 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1533
1534 if (detName == "HLT") {
1535 if (!gROOT->GetClass("AliLevel3")) {
4aa41877 1536 gSystem->Load("libAliHLTSrc.so");
1537 gSystem->Load("libAliHLTMisc.so");
1538 gSystem->Load("libAliHLTHough.so");
1539 gSystem->Load("libAliHLTComp.so");
b8cd5251 1540 }
1541 }
1542
1543 AliReconstructor* reconstructor = NULL;
1544 // first check if a plugin is defined for the reconstructor
1545 TPluginHandler* pluginHandler =
1546 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1547 // if not, add a plugin for it
1548 if (!pluginHandler) {
b8cd5251 1549 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1550 TString libs = gSystem->GetLibraries();
1551 if (libs.Contains("lib" + detName + "base.so") ||
1552 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1553 pluginManager->AddHandler("AliReconstructor", detName,
1554 recName, detName + "rec", recName + "()");
1555 } else {
1556 pluginManager->AddHandler("AliReconstructor", detName,
1557 recName, detName, recName + "()");
c757bafd 1558 }
b8cd5251 1559 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1560 }
1561 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1562 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1563 }
b8cd5251 1564 if (reconstructor) {
1565 TObject* obj = fOptions.FindObject(detName.Data());
1566 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1567 reconstructor->Init(fRunLoader);
b8cd5251 1568 fReconstructor[iDet] = reconstructor;
1569 }
1570
f08fc9f5 1571 // get or create the loader
1572 if (detName != "HLT") {
1573 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1574 if (!fLoader[iDet]) {
1575 AliConfig::Instance()
1576 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1577 detName, detName);
1578 // first check if a plugin is defined for the loader
bb0901a4 1579 pluginHandler =
f08fc9f5 1580 pluginManager->FindHandler("AliLoader", detName);
1581 // if not, add a plugin for it
1582 if (!pluginHandler) {
1583 TString loaderName = "Ali" + detName + "Loader";
1584 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1585 pluginManager->AddHandler("AliLoader", detName,
1586 loaderName, detName + "base",
1587 loaderName + "(const char*, TFolder*)");
1588 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1589 }
1590 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1591 fLoader[iDet] =
1592 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1593 fRunLoader->GetEventFolder());
1594 }
1595 if (!fLoader[iDet]) { // use default loader
1596 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1597 }
1598 if (!fLoader[iDet]) {
1599 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1600 if (fStopOnError) return NULL;
f08fc9f5 1601 } else {
1602 fRunLoader->AddLoader(fLoader[iDet]);
1603 fRunLoader->CdGAFile();
1604 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1605 fRunLoader->Write(0, TObject::kOverwrite);
1606 }
1607 }
1608 }
1609
b8cd5251 1610 return reconstructor;
c757bafd 1611}
1612
2257f27e 1613//_____________________________________________________________________________
1614Bool_t AliReconstruction::CreateVertexer()
1615{
1616// create the vertexer
1617
b8cd5251 1618 fVertexer = NULL;
1619 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1620 if (itsReconstructor) {
b8cd5251 1621 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1622 }
b8cd5251 1623 if (!fVertexer) {
815c2b38 1624 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1625 if (fStopOnError) return kFALSE;
1626 }
1627
1628 return kTRUE;
1629}
1630
24f7a148 1631//_____________________________________________________________________________
b8cd5251 1632Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1633{
f08fc9f5 1634// create the trackers
24f7a148 1635
b8cd5251 1636 TString detStr = detectors;
1637 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1638 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1639 AliReconstructor* reconstructor = GetReconstructor(iDet);
1640 if (!reconstructor) continue;
1641 TString detName = fgkDetectorName[iDet];
1f46a9ae 1642 if (detName == "HLT") {
1643 fRunHLTTracking = kTRUE;
1644 continue;
1645 }
e66fbafb 1646 if (detName == "MUON") {
1647 fRunMuonTracking = kTRUE;
1648 continue;
1649 }
1650
f08fc9f5 1651
1652 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1653 if (!fTracker[iDet] && (iDet < 7)) {
1654 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1655 if (fStopOnError) return kFALSE;
1656 }
1657 }
1658
24f7a148 1659 return kTRUE;
1660}
1661
e583c30d 1662//_____________________________________________________________________________
b26c3770 1663void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1664{
1665// delete trackers and the run loader and close and delete the file
1666
b8cd5251 1667 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1668 delete fReconstructor[iDet];
1669 fReconstructor[iDet] = NULL;
1670 fLoader[iDet] = NULL;
1671 delete fTracker[iDet];
1672 fTracker[iDet] = NULL;
1673 }
1674 delete fVertexer;
1675 fVertexer = NULL;
9178838a 1676 delete fDiamondProfile;
1677 fDiamondProfile = NULL;
e583c30d 1678
1679 delete fRunLoader;
1680 fRunLoader = NULL;
b649205a 1681 delete fRawReader;
1682 fRawReader = NULL;
e583c30d 1683
1684 if (file) {
1685 file->Close();
1686 delete file;
1687 }
b26c3770 1688
1689 if (fileOld) {
1690 fileOld->Close();
1691 delete fileOld;
1692 gSystem->Unlink("AliESDs.old.root");
1693 }
e583c30d 1694}
24f7a148 1695
1696
1697//_____________________________________________________________________________
1698Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1699{
1700// read the ESD event from a file
1701
1702 if (!esd) return kFALSE;
1703 char fileName[256];
1704 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1705 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1706 if (gSystem->AccessPathName(fileName)) return kFALSE;
1707
f3a97c86 1708 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1709 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1710 TFile* file = TFile::Open(fileName);
1711 if (!file || !file->IsOpen()) {
815c2b38 1712 AliError(Form("opening %s failed", fileName));
24f7a148 1713 delete file;
1714 return kFALSE;
1715 }
1716
1717 gROOT->cd();
1718 delete esd;
1719 esd = (AliESD*) file->Get("ESD");
1720 file->Close();
1721 delete file;
1722 return kTRUE;
1723}
1724
1725//_____________________________________________________________________________
1726void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1727{
1728// write the ESD event to a file
1729
1730 if (!esd) return;
1731 char fileName[256];
1732 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1733 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1734
815c2b38 1735 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1736 TFile* file = TFile::Open(fileName, "recreate");
1737 if (!file || !file->IsOpen()) {
815c2b38 1738 AliError(Form("opening %s failed", fileName));
24f7a148 1739 } else {
1740 esd->Write("ESD");
1741 file->Close();
1742 }
1743 delete file;
1744}
f3a97c86 1745
1746
1747
1748
1749//_____________________________________________________________________________
1750void AliReconstruction::CreateTag(TFile* file)
1751{
2bdb9d38 1752 /////////////
1753 //muon code//
1754 ////////////
56982dd1 1755 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1756 //Variables
56982dd1 1757 Double_t fX,fY,fZ ;
1758 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1759 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1760 Int_t fCharge;
1761 TLorentzVector fEPvector;
1762
1763 Float_t fZVertexCut = 10.0;
1764 Float_t fRhoVertexCut = 2.0;
1765
1766 Float_t fLowPtCut = 1.0;
1767 Float_t fHighPtCut = 3.0;
1768 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1769 ////////////
1770
1771 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1772
089bf903 1773 // Creates the tags for all the events in a given ESD file
4302e20f 1774 Int_t ntrack;
089bf903 1775 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1776 Int_t nPos, nNeg, nNeutr;
1777 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1778 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1779 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1780 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1781 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1782 Int_t fVertexflag;
1387d165 1783 Int_t iRunNumber = 0;
547d75a6 1784 TString fVertexName("default");
4302e20f 1785
f3a97c86 1786 AliRunTag *tag = new AliRunTag();
f3a97c86 1787 AliEventTag *evTag = new AliEventTag();
1788 TTree ttag("T","A Tree with event tags");
2bdb9d38 1789 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1790 btag->SetCompressionLevel(9);
2bdb9d38 1791
f3a97c86 1792 AliInfo(Form("Creating the tags......."));
1793
1794 if (!file || !file->IsOpen()) {
1795 AliError(Form("opening failed"));
1796 delete file;
1797 return ;
2bdb9d38 1798 }
d71aa170 1799 Int_t lastEvent = 0;
f3a97c86 1800 TTree *t = (TTree*) file->Get("esdTree");
1801 TBranch * b = t->GetBranch("ESD");
1802 AliESD *esd = 0;
1803 b->SetAddress(&esd);
1387d165 1804
d71aa170 1805 b->GetEntry(fFirstEvent);
1387d165 1806 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1807
bb0901a4 1808 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1809 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1810 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 1811 ntrack = 0;
1812 nPos = 0;
1813 nNeg = 0;
1814 nNeutr =0;
1815 nK0s = 0;
1816 nNeutrons = 0;
1817 nPi0s = 0;
1818 nGamas = 0;
1819 nProtons = 0;
1820 nKaons = 0;
1821 nPions = 0;
1822 nMuons = 0;
1823 nElectrons = 0;
1824 nCh1GeV = 0;
1825 nCh3GeV = 0;
1826 nCh10GeV = 0;
1827 nMu1GeV = 0;
1828 nMu3GeV = 0;
1829 nMu10GeV = 0;
1830 nEl1GeV = 0;
1831 nEl3GeV = 0;
1832 nEl10GeV = 0;
1833 maxPt = .0;
1834 meanPt = .0;
1835 totalP = .0;
547d75a6 1836 fVertexflag = 0;
d1a50cb5 1837
2bdb9d38 1838 b->GetEntry(iEventNumber);
1387d165 1839 iRunNumber = esd->GetRunNumber();
1840 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1841 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1842 if (!vertexIn) AliError("ESD has not defined vertex.");
1843 if (vertexIn) fVertexName = vertexIn->GetName();
1844 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1845 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1846 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1847 UInt_t status = esdTrack->GetStatus();
f3a97c86 1848
2bdb9d38 1849 //select only tracks with ITS refit
1850 if ((status&AliESDtrack::kITSrefit)==0) continue;
1851 //select only tracks with TPC refit
1852 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1853
2bdb9d38 1854 //select only tracks with the "combined PID"
1855 if ((status&AliESDtrack::kESDpid)==0) continue;
1856 Double_t p[3];
1857 esdTrack->GetPxPyPz(p);
1858 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1859 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1860 totalP += momentum;
1861 meanPt += fPt;
1862 if(fPt > maxPt) maxPt = fPt;
4302e20f 1863
2bdb9d38 1864 if(esdTrack->GetSign() > 0) {
1865 nPos++;
56982dd1 1866 if(fPt > fLowPtCut) nCh1GeV++;
1867 if(fPt > fHighPtCut) nCh3GeV++;
1868 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1869 }
1870 if(esdTrack->GetSign() < 0) {
1871 nNeg++;
56982dd1 1872 if(fPt > fLowPtCut) nCh1GeV++;
1873 if(fPt > fHighPtCut) nCh3GeV++;
1874 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1875 }
1876 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1877
2bdb9d38 1878 //PID
1879 Double_t prob[5];
1880 esdTrack->GetESDpid(prob);
4302e20f 1881
2bdb9d38 1882 Double_t rcc = 0.0;
1883 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1884 if(rcc == 0.0) continue;
1885 //Bayes' formula
1886 Double_t w[5];
1887 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1888
2bdb9d38 1889 //protons
1890 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1891 //kaons
1892 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1893 //pions
1894 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1895 //electrons
1896 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1897 nElectrons++;
56982dd1 1898 if(fPt > fLowPtCut) nEl1GeV++;
1899 if(fPt > fHighPtCut) nEl3GeV++;
1900 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1901 }
1902 ntrack++;
1903 }//track loop
1904
1905 /////////////
1906 //muon code//
1907 ////////////
1908 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1909 // loop over all reconstructed tracks (also first track of combination)
1910 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1911 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1912 if (muonTrack == 0x0) continue;
4302e20f 1913
2bdb9d38 1914 // Coordinates at vertex
56982dd1 1915 fZ = muonTrack->GetZ();
1916 fY = muonTrack->GetBendingCoor();
1917 fX = muonTrack->GetNonBendingCoor();
4302e20f 1918
56982dd1 1919 fThetaX = muonTrack->GetThetaX();
1920 fThetaY = muonTrack->GetThetaY();
4302e20f 1921
56982dd1 1922 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1923 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1924 fPxRec = fPzRec * TMath::Tan(fThetaX);
1925 fPyRec = fPzRec * TMath::Tan(fThetaY);
1926 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1927
2bdb9d38 1928 //ChiSquare of the track if needed
56982dd1 1929 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1930 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1931 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1932
2bdb9d38 1933 // total number of muons inside a vertex cut
56982dd1 1934 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1935 nMuons++;
56982dd1 1936 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1937 nMu1GeV++;
56982dd1 1938 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1939 nMu3GeV++;
56982dd1 1940 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1941 nMu10GeV++;
1942 }
1943 }
1944 }
1945 }
1946 }//muon track loop
1947
1948 // Fill the event tags
1949 if(ntrack != 0)
1950 meanPt = meanPt/ntrack;
1951
1952 evTag->SetEventId(iEventNumber+1);
547d75a6 1953 if (vertexIn) {
1954 evTag->SetVertexX(vertexIn->GetXv());
1955 evTag->SetVertexY(vertexIn->GetYv());
1956 evTag->SetVertexZ(vertexIn->GetZv());
1957 evTag->SetVertexZError(vertexIn->GetZRes());
1958 }
d1a50cb5 1959 evTag->SetVertexFlag(fVertexflag);
1960
2bdb9d38 1961 evTag->SetT0VertexZ(esd->GetT0zVertex());
1962
8bd8ac26 1963 evTag->SetTriggerMask(esd->GetTriggerMask());
1964 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1965
32a5cab4 1966 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1967 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1968 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1969 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1970 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1971 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1972
1973
1974 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1975 evTag->SetNumOfPosTracks(nPos);
1976 evTag->SetNumOfNegTracks(nNeg);
1977 evTag->SetNumOfNeutrTracks(nNeutr);
1978
1979 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1980 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1981 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1982 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1983
1984 evTag->SetNumOfProtons(nProtons);
1985 evTag->SetNumOfKaons(nKaons);
1986 evTag->SetNumOfPions(nPions);
1987 evTag->SetNumOfMuons(nMuons);
1988 evTag->SetNumOfElectrons(nElectrons);
1989 evTag->SetNumOfPhotons(nGamas);
1990 evTag->SetNumOfPi0s(nPi0s);
1991 evTag->SetNumOfNeutrons(nNeutrons);
1992 evTag->SetNumOfKaon0s(nK0s);
1993
1994 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1995 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1996 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1997 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1998 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1999 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
2000 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
2001 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
2002 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
2003
85c60a8e 2004 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
2005 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 2006
2007 evTag->SetTotalMomentum(totalP);
2008 evTag->SetMeanPt(meanPt);
2009 evTag->SetMaxPt(maxPt);
2010
1387d165 2011 tag->SetRunId(iInitRunNumber);
2bdb9d38 2012 tag->AddEventTag(*evTag);
2013 }
bb0901a4 2014 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 2015 else lastEvent = fLastEvent;
f3a97c86 2016
2017 ttag.Fill();
2018 tag->Clear();
2019
2020 char fileName[256];
2021 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 2022 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 2023 AliInfo(Form("writing tags to file %s", fileName));
2024 AliDebug(1, Form("writing tags to file %s", fileName));
2025
2026 TFile* ftag = TFile::Open(fileName, "recreate");
2027 ftag->cd();
2028 ttag.Write();
2029 ftag->Close();
2030 file->cd();
2031 delete tag;
f3a97c86 2032 delete evTag;
2033}
2034
a7807689 2035//_____________________________________________________________________________
2036void AliReconstruction::CreateAOD(TFile* esdFile)
2037{
2038 // do nothing for now
2039
2040 return;
2041}
2042
2043
98937d93 2044void AliReconstruction::WriteAlignmentData(AliESD* esd)
2045{
2046 // Write space-points which are then used in the alignment procedures
2047 // For the moment only ITS, TRD and TPC
2048
2049 // Load TOF clusters
d528ee75 2050 if (fTracker[3]){
2051 fLoader[3]->LoadRecPoints("read");
2052 TTree* tree = fLoader[3]->TreeR();
2053 if (!tree) {
2054 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2055 return;
2056 }
2057 fTracker[3]->LoadClusters(tree);
98937d93 2058 }
98937d93 2059 Int_t ntracks = esd->GetNumberOfTracks();
2060 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2061 {
2062 AliESDtrack *track = esd->GetTrack(itrack);
2063 Int_t nsp = 0;
ef7253ac 2064 Int_t idx[200];
98937d93 2065 for (Int_t iDet = 3; iDet >= 0; iDet--)
2066 nsp += track->GetNcls(iDet);
2067 if (nsp) {
2068 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2069 track->SetTrackPointArray(sp);
2070 Int_t isptrack = 0;
2071 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2072 AliTracker *tracker = fTracker[iDet];
2073 if (!tracker) continue;
2074 Int_t nspdet = track->GetNcls(iDet);
98937d93 2075 if (nspdet <= 0) continue;
2076 track->GetClusters(iDet,idx);
2077 AliTrackPoint p;
2078 Int_t isp = 0;
2079 Int_t isp2 = 0;
2080 while (isp < nspdet) {
2081 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 2082 const Int_t kNTPCmax = 159;
2083 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2084 if (!isvalid) continue;
2085 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2086 }
98937d93 2087 }
2088 }
2089 }
d528ee75 2090 if (fTracker[3]){
2091 fTracker[3]->UnloadClusters();
2092 fLoader[3]->UnloadRecPoints();
2093 }
98937d93 2094}
2e3550da 2095
2096//_____________________________________________________________________________
2097void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2098{
2099 // The method reads the raw-data error log
2100 // accumulated within the rawReader.
2101 // It extracts the raw-data errors related to
2102 // the current event and stores them into
2103 // a TClonesArray inside the esd object.
2104
2105 if (!fRawReader) return;
2106
2107 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2108
2109 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2110 if (!log) continue;
2111 if (iEvent != log->GetEventNumber()) continue;
2112
2113 esd->AddRawDataErrorLog(log);
2114 }
2115
2116}