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