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