]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliReconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Rename AliTRDCalGlobals to AliTRDCalFEE and remove fADCClockphase and fPretriggerConf
[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
25be1e5c 154#include "AliGeomManager.h"
98937d93 155#include "AliTrackPointArray.h"
b0314964 156#include "AliCDBManager.h"
6bae477a 157#include "AliCDBEntry.h"
158#include "AliAlignObj.h"
f3a97c86 159
b647652d 160#include "AliCentralTrigger.h"
161#include "AliCTPRawStream.h"
162
f29f1726 163#include "AliAODEvent.h"
164#include "AliAODHeader.h"
165#include "AliAODTrack.h"
166#include "AliAODVertex.h"
167#include "AliAODCluster.h"
168
169
596a855f 170ClassImp(AliReconstruction)
171
172
c757bafd 173//_____________________________________________________________________________
b384f8a4 174const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
c757bafd 175
596a855f 176//_____________________________________________________________________________
024cf675 177AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 178 const char* name, const char* title) :
179 TNamed(name, title),
180
c84a5e9e 181 fUniformField(kTRUE),
2257f27e 182 fRunVertexFinder(kTRUE),
1f46a9ae 183 fRunHLTTracking(kFALSE),
e66fbafb 184 fRunMuonTracking(kFALSE),
1d99986f 185 fStopOnError(kFALSE),
186 fWriteAlignmentData(kFALSE),
187 fWriteESDfriend(kFALSE),
a7807689 188 fWriteAOD(kFALSE),
b647652d 189 fFillTriggerESD(kTRUE),
1d99986f 190
191 fRunLocalReconstruction("ALL"),
b8cd5251 192 fRunTracking("ALL"),
e583c30d 193 fFillESD("ALL"),
194 fGAliceFileName(gAliceFilename),
b649205a 195 fInput(""),
35042093 196 fEquipIdMap(""),
b26c3770 197 fFirstEvent(0),
198 fLastEvent(-1),
973388c2 199 fNumberOfEventsPerFile(1),
24f7a148 200 fCheckPointLevel(0),
b8cd5251 201 fOptions(),
6bae477a 202 fLoadAlignFromCDB(kTRUE),
203 fLoadAlignData("ALL"),
46698ae4 204 fESDPar(""),
e583c30d 205
206 fRunLoader(NULL),
b649205a 207 fRawReader(NULL),
b8cd5251 208
98937d93 209 fVertexer(NULL),
9178838a 210 fDiamondProfile(NULL),
98937d93 211
6bae477a 212 fAlignObjArray(NULL),
ec92bee0 213 fCDBUri(cdbUri),
214 fSpecCDBUri()
596a855f 215{
216// create reconstruction object with default parameters
b8cd5251 217
218 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
219 fReconstructor[iDet] = NULL;
220 fLoader[iDet] = NULL;
221 fTracker[iDet] = NULL;
222 }
e47c4c2e 223 AliPID pid;
596a855f 224}
225
226//_____________________________________________________________________________
227AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 228 TNamed(rec),
229
c84a5e9e 230 fUniformField(rec.fUniformField),
2257f27e 231 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 232 fRunHLTTracking(rec.fRunHLTTracking),
e66fbafb 233 fRunMuonTracking(rec.fRunMuonTracking),
1d99986f 234 fStopOnError(rec.fStopOnError),
235 fWriteAlignmentData(rec.fWriteAlignmentData),
236 fWriteESDfriend(rec.fWriteESDfriend),
a7807689 237 fWriteAOD(rec.fWriteAOD),
b647652d 238 fFillTriggerESD(rec.fFillTriggerESD),
1d99986f 239
240 fRunLocalReconstruction(rec.fRunLocalReconstruction),
e583c30d 241 fRunTracking(rec.fRunTracking),
242 fFillESD(rec.fFillESD),
243 fGAliceFileName(rec.fGAliceFileName),
b649205a 244 fInput(rec.fInput),
35042093 245 fEquipIdMap(rec.fEquipIdMap),
b26c3770 246 fFirstEvent(rec.fFirstEvent),
247 fLastEvent(rec.fLastEvent),
973388c2 248 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
24f7a148 249 fCheckPointLevel(0),
b8cd5251 250 fOptions(),
6bae477a 251 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
252 fLoadAlignData(rec.fLoadAlignData),
46698ae4 253 fESDPar(rec.fESDPar),
e583c30d 254
255 fRunLoader(NULL),
b649205a 256 fRawReader(NULL),
b8cd5251 257
98937d93 258 fVertexer(NULL),
9178838a 259 fDiamondProfile(NULL),
98937d93 260
6bae477a 261 fAlignObjArray(rec.fAlignObjArray),
ec92bee0 262 fCDBUri(rec.fCDBUri),
263 fSpecCDBUri()
596a855f 264{
265// copy constructor
266
ec92bee0 267 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
efd2085e 268 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
269 }
b8cd5251 270 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
271 fReconstructor[iDet] = NULL;
272 fLoader[iDet] = NULL;
273 fTracker[iDet] = NULL;
274 }
ec92bee0 275 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
276 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
277 }
596a855f 278}
279
280//_____________________________________________________________________________
281AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
282{
283// assignment operator
284
285 this->~AliReconstruction();
286 new(this) AliReconstruction(rec);
287 return *this;
288}
289
290//_____________________________________________________________________________
291AliReconstruction::~AliReconstruction()
292{
293// clean up
294
e583c30d 295 CleanUp();
efd2085e 296 fOptions.Delete();
ec92bee0 297 fSpecCDBUri.Delete();
596a855f 298}
299
024cf675 300//_____________________________________________________________________________
301void AliReconstruction::InitCDBStorage()
302{
303// activate a default CDB storage
304// First check if we have any CDB storage set, because it is used
305// to retrieve the calibration and alignment constants
306
307 AliCDBManager* man = AliCDBManager::Instance();
ec92bee0 308 if (man->IsDefaultStorageSet())
024cf675 309 {
310 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 311 AliWarning("Default CDB storage has been already set !");
312 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
024cf675 313 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 314 fCDBUri = "";
315 }
316 else {
b8ec52f6 317 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
318 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
319 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 320 man->SetDefaultStorage(fCDBUri);
321 }
322
323 // Now activate the detector specific CDB storage locations
c3a7b59a 324 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
325 TObject* obj = fSpecCDBUri[i];
326 if (!obj) continue;
b8ec52f6 327 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
328 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
329 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
c3a7b59a 330 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
ec92bee0 331 }
727d0766 332 man->Print();
024cf675 333}
334
335//_____________________________________________________________________________
336void AliReconstruction::SetDefaultStorage(const char* uri) {
ec92bee0 337// Store the desired default CDB storage location
338// Activate it later within the Run() method
024cf675 339
ec92bee0 340 fCDBUri = uri;
024cf675 341
342}
343
344//_____________________________________________________________________________
c3a7b59a 345void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
ec92bee0 346// Store a detector-specific CDB storage location
347// Activate it later within the Run() method
024cf675 348
c3a7b59a 349 AliCDBPath aPath(calibType);
350 if(!aPath.IsValid()){
351 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
352 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
353 if(!strcmp(calibType, fgkDetectorName[iDet])) {
354 aPath.SetPath(Form("%s/*", calibType));
355 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
356 break;
357 }
358 }
359 if(!aPath.IsValid()){
360 AliError(Form("Not a valid path or detector: %s", calibType));
361 return;
362 }
363 }
364
365 // check that calibType refers to a "valid" detector name
366 Bool_t isDetector = kFALSE;
367 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
368 TString detName = fgkDetectorName[iDet];
369 if(aPath.GetLevel0() == detName) {
370 isDetector = kTRUE;
371 break;
372 }
373 }
374
375 if(!isDetector) {
376 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
377 return;
378 }
379
380 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
ec92bee0 381 if (obj) fSpecCDBUri.Remove(obj);
c3a7b59a 382 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
024cf675 383
384}
385
46698ae4 386
387
388
6bae477a 389//_____________________________________________________________________________
390Bool_t AliReconstruction::SetRunNumber()
391{
392 // The method is called in Run() in order
393 // to set a correct run number.
394 // In case of raw data reconstruction the
395 // run number is taken from the raw data header
396
397 if(AliCDBManager::Instance()->GetRun() < 0) {
398 if (!fRunLoader) {
399 AliError("No run loader is found !");
400 return kFALSE;
401 }
402 // read run number from gAlice
ec92bee0 403 if(fRunLoader->GetAliRun())
404 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
405 else {
406 if(fRawReader) {
407 if(fRawReader->NextEvent()) {
408 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
409 fRawReader->RewindEvents();
410 }
411 else {
412 AliError("No raw-data events found !");
413 return kFALSE;
414 }
415 }
416 else {
417 AliError("Neither gAlice nor RawReader objects are found !");
418 return kFALSE;
419 }
420 }
421 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
6bae477a 422 }
423 return kTRUE;
424}
425
6bae477a 426//_____________________________________________________________________________
427Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
428{
429 // Read the alignment objects from CDB.
430 // Each detector is supposed to have the
431 // alignment objects in DET/Align/Data CDB path.
432 // All the detector objects are then collected,
433 // sorted by geometry level (starting from ALIC) and
434 // then applied to the TGeo geometry.
435 // Finally an overlaps check is performed.
436
437 // Load alignment data from CDB and fill fAlignObjArray
438 if(fLoadAlignFromCDB){
6bae477a 439
25be1e5c 440 TString detStr = detectors;
441 TString loadAlObjsListOfDets = "";
442
443 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
444 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
445 loadAlObjsListOfDets += fgkDetectorName[iDet];
446 loadAlObjsListOfDets += " ";
447 } // end loop over detectors
448 (AliGeomManager::Instance())->ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
449 }else{
450 // Check if the array with alignment objects was
451 // provided by the user. If yes, apply the objects
452 // to the present TGeo geometry
453 if (fAlignObjArray) {
454 if (gGeoManager && gGeoManager->IsClosed()) {
455 if ((AliGeomManager::Instance())->ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
456 AliError("The misalignment of one or more volumes failed!"
457 "Compare the list of simulated detectors and the list of detector alignment data!");
458 return kFALSE;
459 }
460 }
461 else {
462 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
6bae477a 463 return kFALSE;
464 }
465 }
6bae477a 466 }
25be1e5c 467
8e245d15 468 delete fAlignObjArray; fAlignObjArray=0;
a03b0371 469
470 // Update the TGeoPhysicalNodes
471 gGeoManager->RefreshPhysicalNodes();
472
6bae477a 473 return kTRUE;
474}
596a855f 475
476//_____________________________________________________________________________
477void AliReconstruction::SetGAliceFile(const char* fileName)
478{
479// set the name of the galice file
480
481 fGAliceFileName = fileName;
482}
483
efd2085e 484//_____________________________________________________________________________
485void AliReconstruction::SetOption(const char* detector, const char* option)
486{
487// set options for the reconstruction of a detector
488
489 TObject* obj = fOptions.FindObject(detector);
490 if (obj) fOptions.Remove(obj);
491 fOptions.Add(new TNamed(detector, option));
492}
493
596a855f 494
495//_____________________________________________________________________________
4a33489c 496Bool_t AliReconstruction::Run(const char* input)
596a855f 497{
498// run the reconstruction
499
b649205a 500 // set the input
501 if (!input) input = fInput.Data();
502 TString fileName(input);
503 if (fileName.EndsWith("/")) {
504 fRawReader = new AliRawReaderFile(fileName);
505 } else if (fileName.EndsWith(".root")) {
506 fRawReader = new AliRawReaderRoot(fileName);
507 } else if (!fileName.IsNull()) {
508 fRawReader = new AliRawReaderDate(fileName);
509 fRawReader->SelectEvents(7);
510 }
35042093 511 if (!fEquipIdMap.IsNull() && fRawReader)
512 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
513
b649205a 514
f08fc9f5 515 // get the run loader
516 if (!InitRunLoader()) return kFALSE;
596a855f 517
ec92bee0 518 // Initialize the CDB storage
519 InitCDBStorage();
520
6bae477a 521 // Set run number in CDBManager (if it is not already set by the user)
522 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
523
524 // Import ideal TGeo geometry and apply misalignment
525 if (!gGeoManager) {
526 TString geom(gSystem->DirName(fGAliceFileName));
527 geom += "/geometry.root";
528 TGeoManager::Import(geom.Data());
529 if (!gGeoManager) if (fStopOnError) return kFALSE;
530 }
8e245d15 531
532 AliCDBManager* man = AliCDBManager::Instance();
6bae477a 533 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
534
596a855f 535 // local reconstruction
59697224 536 if (!fRunLocalReconstruction.IsNull()) {
537 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 538 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 539 }
540 }
b26c3770 541// if (!fRunVertexFinder && fRunTracking.IsNull() &&
542// fFillESD.IsNull()) return kTRUE;
2257f27e 543
544 // get vertexer
545 if (fRunVertexFinder && !CreateVertexer()) {
546 if (fStopOnError) {
547 CleanUp();
548 return kFALSE;
549 }
550 }
596a855f 551
f08fc9f5 552 // get trackers
b8cd5251 553 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 554 if (fStopOnError) {
555 CleanUp();
556 return kFALSE;
557 }
596a855f 558 }
24f7a148 559
9db3a215 560
561 TStopwatch stopwatch;
562 stopwatch.Start();
563
b26c3770 564 // get the possibly already existing ESD file and tree
46698ae4 565 AliESD* esd = new AliESD(); AliESD* hltesd = new AliESD();
b26c3770 566 TFile* fileOld = NULL;
1f46a9ae 567 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 568 if (!gSystem->AccessPathName("AliESDs.root")){
569 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
570 fileOld = TFile::Open("AliESDs.old.root");
571 if (fileOld && fileOld->IsOpen()) {
572 treeOld = (TTree*) fileOld->Get("esdTree");
46698ae4 573 if (treeOld)esd->ReadFromTree(treeOld);
1f46a9ae 574 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
46698ae4 575 if (hlttreeOld) hltesd->ReadFromTree(hlttreeOld);
b26c3770 576 }
577 }
578
36711aa4 579 // create the ESD output file and tree
596a855f 580 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
46698ae4 581 file->SetCompressionLevel(2);
596a855f 582 if (!file->IsOpen()) {
815c2b38 583 AliError("opening AliESDs.root failed");
b26c3770 584 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 585 }
46698ae4 586
36711aa4 587 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
46698ae4 588 esd = new AliESD();
589 esd->CreateStdContent();
590 esd->WriteToTree(tree);
591
1f46a9ae 592 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
46698ae4 593 hltesd = new AliESD();
594 hltesd->CreateStdContent();
595 hltesd->WriteToTree(hlttree);
596
597 /* CKB Why?
1f46a9ae 598 delete esd; delete hltesd;
599 esd = NULL; hltesd = NULL;
46698ae4 600 */
500d54ab 601 // create the branch with ESD additions
46698ae4 602 AliESDfriend *esdf = 0;
1d99986f 603 if (fWriteESDfriend) {
46698ae4 604 esdf = new AliESDfriend();
605 TBranch *br=tree->Branch("ESDfriend.","AliESDfriend", &esdf);
606 br->SetFile("AliESDfriends.root");
607 esd->AddObject(esdf);
1d99986f 608 }
46698ae4 609
17c86e90 610 // Get the diamond profile from OCDB
611 AliCDBEntry* entry = AliCDBManager::Instance()
612 ->Get("GRP/Calib/MeanVertex");
613
614 if(entry) {
615 fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
616 } else {
617 AliError("No diamond profile found in OCDB!");
618 }
619
20e5681c 620 AliVertexerTracks tVertexer(AliTracker::GetBz());
9178838a 621 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
c5e3e5d1 622
596a855f 623 // loop over events
b649205a 624 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 625
596a855f 626 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 627 if (fRawReader) fRawReader->NextEvent();
4a33489c 628 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
b26c3770 629 // copy old ESD to the new one
630 if (treeOld) {
46698ae4 631 esd->ReadFromTree(treeOld);
b26c3770 632 treeOld->GetEntry(iEvent);
633 }
634 tree->Fill();
1f46a9ae 635 if (hlttreeOld) {
46698ae4 636 esd->ReadFromTree(hlttreeOld);
1f46a9ae 637 hlttreeOld->GetEntry(iEvent);
638 }
639 hlttree->Fill();
b26c3770 640 continue;
641 }
46698ae4 642
815c2b38 643 AliInfo(Form("processing event %d", iEvent));
596a855f 644 fRunLoader->GetEvent(iEvent);
24f7a148 645
bb0901a4 646 char aFileName[256];
647 sprintf(aFileName, "ESD_%d.%d_final.root",
f08fc9f5 648 fRunLoader->GetHeader()->GetRun(),
649 fRunLoader->GetHeader()->GetEventNrInRun());
bb0901a4 650 if (!gSystem->AccessPathName(aFileName)) continue;
24f7a148 651
b26c3770 652 // local reconstruction
653 if (!fRunLocalReconstruction.IsNull()) {
654 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
655 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
656 }
657 }
658
46698ae4 659
f08fc9f5 660 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 661 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
31fd97b2 662 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
663 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
46698ae4 664
d6ee376f 665 // Set magnetic field from the tracker
666 esd->SetMagneticField(AliTracker::GetBz());
667 hltesd->SetMagneticField(AliTracker::GetBz());
596a855f 668
46698ae4 669
670
2e3550da 671 // Fill raw-data error log into the ESD
672 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
673
2257f27e 674 // vertex finder
675 if (fRunVertexFinder) {
676 if (!ReadESD(esd, "vertex")) {
677 if (!RunVertexFinder(esd)) {
b26c3770 678 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 679 }
680 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
681 }
682 }
683
1f46a9ae 684 // HLT tracking
685 if (!fRunTracking.IsNull()) {
686 if (fRunHLTTracking) {
687 hltesd->SetVertex(esd->GetVertex());
688 if (!RunHLTTracking(hltesd)) {
689 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
690 }
691 }
692 }
693
e66fbafb 694 // Muon tracking
b8cd5251 695 if (!fRunTracking.IsNull()) {
e66fbafb 696 if (fRunMuonTracking) {
697 if (!RunMuonTracking()) {
b26c3770 698 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 699 }
596a855f 700 }
701 }
702
e66fbafb 703 // barrel tracking
704 if (!fRunTracking.IsNull()) {
21c573b7 705 if (!ReadESD(esd, "tracking")) {
706 if (!RunTracking(esd)) {
707 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
e66fbafb 708 }
21c573b7 709 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
e66fbafb 710 }
711 }
21c573b7 712
596a855f 713 // fill ESD
714 if (!fFillESD.IsNull()) {
715 if (!FillESD(esd, fFillESD)) {
b26c3770 716 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 717 }
718 }
001397cd 719 // fill Event header information from the RawEventHeader
720 if (fRawReader){FillRawEventHeaderESD(esd);}
596a855f 721
722 // combined PID
723 AliESDpid::MakePID(esd);
24f7a148 724 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 725
b647652d 726 if (fFillTriggerESD) {
727 if (!ReadESD(esd, "trigger")) {
728 if (!FillTriggerESD(esd)) {
729 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
730 }
731 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
732 }
733 }
734
a03cd2e0 735 //Try to improve the reconstructed primary vertex position using the tracks
736 AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
17c86e90 737 if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
738
a03cd2e0 739 if (pvtx)
740 if (pvtx->GetStatus()) {
741 // Store the improved primary vertex
742 esd->SetPrimaryVertex(pvtx);
743 // Propagate the tracks to the DCA to the improved primary vertex
744 Double_t somethingbig = 777.;
745 Double_t bz = esd->GetMagneticField();
746 Int_t nt=esd->GetNumberOfTracks();
747 while (nt--) {
748 AliESDtrack *t = esd->GetTrack(nt);
749 t->RelateToVertex(pvtx, bz, somethingbig);
750 }
751 }
c5e3e5d1 752
5e4ff34d 753 {
754 // V0 finding
755 AliV0vertexer vtxer;
756 vtxer.Tracks2V0vertices(esd);
757
758 // Cascade finding
759 AliCascadeVertexer cvtxer;
760 cvtxer.V0sTracks2CascadeVertices(esd);
761 }
762
596a855f 763 // write ESD
1d99986f 764 if (fWriteESDfriend) {
46698ae4 765 new (esdf) AliESDfriend(); // Reset...
766 esd->GetESDfriend(esdf);
1d99986f 767 }
500d54ab 768 tree->Fill();
769
770 // write HLT ESD
771 hlttree->Fill();
1d99986f 772
f3a97c86 773 if (fCheckPointLevel > 0) WriteESD(esd, "final");
46698ae4 774 esd->Reset();
775 hltesd->Reset();
776 // esdf->Reset();
777 // delete esdf; esdf = 0;
778 }
779
780
781 tree->GetUserInfo()->Add(esd);
782 hlttree->GetUserInfo()->Add(hltesd);
783
784
785
786 if(fESDPar.Contains("ESD.par")){
787 AliInfo("Attaching ESD.par to Tree");
788 TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
789 tree->GetUserInfo()->Add(fn);
596a855f 790 }
791
46698ae4 792
9db3a215 793 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
794 stopwatch.RealTime(),stopwatch.CpuTime()));
795
36711aa4 796 file->cd();
a9c0e6db 797 if (fWriteESDfriend)
798 tree->SetBranchStatus("ESDfriend*",0);
36711aa4 799 tree->Write();
1f46a9ae 800 hlttree->Write();
f3a97c86 801
a7807689 802 if (fWriteAOD) {
f29f1726 803 TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
804 ESDFile2AODFile(file, aodFile);
805 aodFile->Close();
a7807689 806 }
807
46698ae4 808 gROOT->cd();
f3a97c86 809 // Create tags for the events in the ESD tree (the ESD tree is always present)
810 // In case of empty events the tags will contain dummy values
811 CreateTag(file);
46698ae4 812
813 CleanUp(file, fileOld);
596a855f 814
815 return kTRUE;
816}
817
818
819//_____________________________________________________________________________
59697224 820Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 821{
59697224 822// run the local reconstruction
596a855f 823
030b532d 824 TStopwatch stopwatch;
825 stopwatch.Start();
826
8e245d15 827 AliCDBManager* man = AliCDBManager::Instance();
828 Bool_t origCache = man->GetCacheFlag();
829
596a855f 830 TString detStr = detectors;
b8cd5251 831 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
832 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
833 AliReconstructor* reconstructor = GetReconstructor(iDet);
834 if (!reconstructor) continue;
b26c3770 835 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 836
837 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
838 TStopwatch stopwatchDet;
839 stopwatchDet.Start();
8e245d15 840
841 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
842
843 man->SetCacheFlag(kTRUE);
844 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
845 man->GetAll(calibPath); // entries are cached!
846
b8cd5251 847 if (fRawReader) {
848 fRawReader->RewindEvents();
849 reconstructor->Reconstruct(fRunLoader, fRawReader);
850 } else {
851 reconstructor->Reconstruct(fRunLoader);
596a855f 852 }
5f8272e1 853 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
854 fgkDetectorName[iDet],
855 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
8e245d15 856
857 // unload calibration data
858 man->ClearCache();
596a855f 859 }
860
8e245d15 861 man->SetCacheFlag(origCache);
862
596a855f 863 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 864 AliError(Form("the following detectors were not found: %s",
865 detStr.Data()));
596a855f 866 if (fStopOnError) return kFALSE;
867 }
868
5f8272e1 869 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
870 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 871
596a855f 872 return kTRUE;
873}
874
b26c3770 875//_____________________________________________________________________________
876Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
877{
878// run the local reconstruction
879
880 TStopwatch stopwatch;
881 stopwatch.Start();
882
883 TString detStr = detectors;
884 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
885 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
886 AliReconstructor* reconstructor = GetReconstructor(iDet);
887 if (!reconstructor) continue;
888 AliLoader* loader = fLoader[iDet];
889
890 // conversion of digits
891 if (fRawReader && reconstructor->HasDigitConversion()) {
892 AliInfo(Form("converting raw data digits into root objects for %s",
893 fgkDetectorName[iDet]));
894 TStopwatch stopwatchDet;
895 stopwatchDet.Start();
896 loader->LoadDigits("update");
897 loader->CleanDigits();
898 loader->MakeDigitsContainer();
899 TTree* digitsTree = loader->TreeD();
900 reconstructor->ConvertDigits(fRawReader, digitsTree);
901 loader->WriteDigits("OVERWRITE");
902 loader->UnloadDigits();
5f8272e1 903 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
904 fgkDetectorName[iDet],
905 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 906 }
907
908 // local reconstruction
909 if (!reconstructor->HasLocalReconstruction()) continue;
910 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
911 TStopwatch stopwatchDet;
912 stopwatchDet.Start();
913 loader->LoadRecPoints("update");
914 loader->CleanRecPoints();
915 loader->MakeRecPointsContainer();
916 TTree* clustersTree = loader->TreeR();
917 if (fRawReader && !reconstructor->HasDigitConversion()) {
918 reconstructor->Reconstruct(fRawReader, clustersTree);
919 } else {
920 loader->LoadDigits("read");
921 TTree* digitsTree = loader->TreeD();
922 if (!digitsTree) {
923 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
924 if (fStopOnError) return kFALSE;
925 } else {
926 reconstructor->Reconstruct(digitsTree, clustersTree);
927 }
928 loader->UnloadDigits();
929 }
930 loader->WriteRecPoints("OVERWRITE");
931 loader->UnloadRecPoints();
ade5466c 932 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
933 fgkDetectorName[iDet],
934 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 935 }
936
937 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
938 AliError(Form("the following detectors were not found: %s",
939 detStr.Data()));
940 if (fStopOnError) return kFALSE;
941 }
5f8272e1 942
943 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
944 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 945
946 return kTRUE;
947}
948
596a855f 949//_____________________________________________________________________________
2257f27e 950Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 951{
952// run the barrel tracking
953
030b532d 954 TStopwatch stopwatch;
955 stopwatch.Start();
956
2257f27e 957 AliESDVertex* vertex = NULL;
958 Double_t vtxPos[3] = {0, 0, 0};
959 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
960 TArrayF mcVertex(3);
a6b0b91b 961 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
962 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
963 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
964 }
2257f27e 965
b8cd5251 966 if (fVertexer) {
17c86e90 967 if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
815c2b38 968 AliInfo("running the ITS vertex finder");
b26c3770 969 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 970 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 971 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 972 if(!vertex){
815c2b38 973 AliWarning("Vertex not found");
c710f220 974 vertex = new AliESDVertex();
d1a50cb5 975 vertex->SetName("default");
2257f27e 976 }
977 else {
d1a50cb5 978 vertex->SetName("reconstructed");
2257f27e 979 }
980
981 } else {
815c2b38 982 AliInfo("getting the primary vertex from MC");
2257f27e 983 vertex = new AliESDVertex(vtxPos, vtxErr);
984 }
985
986 if (vertex) {
987 vertex->GetXYZ(vtxPos);
988 vertex->GetSigmaXYZ(vtxErr);
989 } else {
815c2b38 990 AliWarning("no vertex reconstructed");
2257f27e 991 vertex = new AliESDVertex(vtxPos, vtxErr);
992 }
993 esd->SetVertex(vertex);
32e449be 994 // if SPD multiplicity has been determined, it is stored in the ESD
25be1e5c 995 AliMultiplicity *mult = fVertexer->GetMultiplicity();
32e449be 996 if(mult)esd->SetMultiplicity(mult);
997
b8cd5251 998 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
999 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1000 }
2257f27e 1001 delete vertex;
1002
5f8272e1 1003 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1004 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 1005
1006 return kTRUE;
1007}
1008
1f46a9ae 1009//_____________________________________________________________________________
1010Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1011{
1012// run the HLT barrel tracking
1013
1014 TStopwatch stopwatch;
1015 stopwatch.Start();
1016
1017 if (!fRunLoader) {
1018 AliError("Missing runLoader!");
1019 return kFALSE;
1020 }
1021
1022 AliInfo("running HLT tracking");
1023
1024 // Get a pointer to the HLT reconstructor
1025 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1026 if (!reconstructor) return kFALSE;
1027
1028 // TPC + ITS
1029 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1030 TString detName = fgkDetectorName[iDet];
1031 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1032 reconstructor->SetOption(detName.Data());
1033 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1034 if (!tracker) {
1035 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1036 if (fStopOnError) return kFALSE;
9dcc06e1 1037 continue;
1f46a9ae 1038 }
1039 Double_t vtxPos[3];
1040 Double_t vtxErr[3]={0.005,0.005,0.010};
1041 const AliESDVertex *vertex = esd->GetVertex();
1042 vertex->GetXYZ(vtxPos);
1043 tracker->SetVertex(vtxPos,vtxErr);
1044 if(iDet != 1) {
1045 fLoader[iDet]->LoadRecPoints("read");
1046 TTree* tree = fLoader[iDet]->TreeR();
1047 if (!tree) {
1048 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1049 return kFALSE;
1050 }
1051 tracker->LoadClusters(tree);
1052 }
1053 if (tracker->Clusters2Tracks(esd) != 0) {
1054 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1055 return kFALSE;
1056 }
1057 if(iDet != 1) {
1058 tracker->UnloadClusters();
1059 }
1060 delete tracker;
1061 }
1062
5f8272e1 1063 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1064 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 1065
1066 return kTRUE;
1067}
1068
e66fbafb 1069//_____________________________________________________________________________
1070Bool_t AliReconstruction::RunMuonTracking()
1071{
1072// run the muon spectrometer tracking
1073
1074 TStopwatch stopwatch;
1075 stopwatch.Start();
1076
1077 if (!fRunLoader) {
1078 AliError("Missing runLoader!");
1079 return kFALSE;
1080 }
1081 Int_t iDet = 7; // for MUON
1082
1083 AliInfo("is running...");
1084
1085 // Get a pointer to the MUON reconstructor
1086 AliReconstructor *reconstructor = GetReconstructor(iDet);
1087 if (!reconstructor) return kFALSE;
1088
1089
1090 TString detName = fgkDetectorName[iDet];
1091 AliDebug(1, Form("%s tracking", detName.Data()));
1092 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1093 if (!tracker) {
1094 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1095 return kFALSE;
1096 }
1097
1098 // create Tracks
1099 fLoader[iDet]->LoadTracks("update");
1100 fLoader[iDet]->CleanTracks();
1101 fLoader[iDet]->MakeTracksContainer();
1102
1103 // read RecPoints
1104 fLoader[iDet]->LoadRecPoints("read");
1105
1106 if (!tracker->Clusters2Tracks(0x0)) {
1107 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1108 return kFALSE;
1109 }
1110 fLoader[iDet]->UnloadRecPoints();
1111
1112 fLoader[iDet]->WriteTracks("OVERWRITE");
1113 fLoader[iDet]->UnloadTracks();
1114
1115 delete tracker;
1116
1117
1118 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1119 stopwatch.RealTime(),stopwatch.CpuTime()));
1120
1121 return kTRUE;
1122}
1123
1124
2257f27e 1125//_____________________________________________________________________________
1126Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1127{
1128// run the barrel tracking
1129
1130 TStopwatch stopwatch;
1131 stopwatch.Start();
24f7a148 1132
815c2b38 1133 AliInfo("running tracking");
596a855f 1134
91b876d1 1135 //Fill the ESD with the T0 info (will be used by the TOF)
1f4331b3 1136 if (fReconstructor[11])
1137 GetReconstructor(11)->FillESD(fRunLoader, esd);
91b876d1 1138
b8cd5251 1139 // pass 1: TPC + ITS inwards
1140 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1141 if (!fTracker[iDet]) continue;
1142 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 1143
b8cd5251 1144 // load clusters
1145 fLoader[iDet]->LoadRecPoints("read");
1146 TTree* tree = fLoader[iDet]->TreeR();
1147 if (!tree) {
1148 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1149 return kFALSE;
1150 }
b8cd5251 1151 fTracker[iDet]->LoadClusters(tree);
1152
1153 // run tracking
1154 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1155 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 1156 return kFALSE;
1157 }
b8cd5251 1158 if (fCheckPointLevel > 1) {
1159 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1160 }
878e1fe1 1161 // preliminary PID in TPC needed by the ITS tracker
1162 if (iDet == 1) {
1163 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 1164 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 1165 AliESDpid::MakePID(esd);
1166 }
b8cd5251 1167 }
596a855f 1168
b8cd5251 1169 // pass 2: ALL backwards
1170 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1171 if (!fTracker[iDet]) continue;
1172 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1173
1174 // load clusters
1175 if (iDet > 1) { // all except ITS, TPC
1176 TTree* tree = NULL;
7b61cd9c 1177 fLoader[iDet]->LoadRecPoints("read");
1178 tree = fLoader[iDet]->TreeR();
b8cd5251 1179 if (!tree) {
1180 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1181 return kFALSE;
1182 }
b8cd5251 1183 fTracker[iDet]->LoadClusters(tree);
1184 }
24f7a148 1185
b8cd5251 1186 // run tracking
1187 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1188 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
49dfd67a 1189 // return kFALSE;
b8cd5251 1190 }
1191 if (fCheckPointLevel > 1) {
1192 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1193 }
24f7a148 1194
b8cd5251 1195 // unload clusters
1196 if (iDet > 2) { // all except ITS, TPC, TRD
1197 fTracker[iDet]->UnloadClusters();
7b61cd9c 1198 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1199 }
8f37df88 1200 // updated PID in TPC needed by the ITS tracker -MI
1201 if (iDet == 1) {
1202 GetReconstructor(1)->FillESD(fRunLoader, esd);
1203 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1204 AliESDpid::MakePID(esd);
1205 }
b8cd5251 1206 }
596a855f 1207
98937d93 1208 // write space-points to the ESD in case alignment data output
1209 // is switched on
1210 if (fWriteAlignmentData)
1211 WriteAlignmentData(esd);
1212
b8cd5251 1213 // pass 3: TRD + TPC + ITS refit inwards
1214 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1215 if (!fTracker[iDet]) continue;
1216 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 1217
b8cd5251 1218 // run tracking
1219 if (fTracker[iDet]->RefitInward(esd) != 0) {
1220 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
49dfd67a 1221 // return kFALSE;
b8cd5251 1222 }
1223 if (fCheckPointLevel > 1) {
1224 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1225 }
596a855f 1226
b8cd5251 1227 // unload clusters
1228 fTracker[iDet]->UnloadClusters();
1229 fLoader[iDet]->UnloadRecPoints();
1230 }
ff8bb5ae 1231 //
1232 // Propagate track to the vertex - if not done by ITS
1233 //
1234 Int_t ntracks = esd->GetNumberOfTracks();
1235 for (Int_t itrack=0; itrack<ntracks; itrack++){
1236 const Double_t kRadius = 3; // beam pipe radius
1237 const Double_t kMaxStep = 5; // max step
1238 const Double_t kMaxD = 123456; // max distance to prim vertex
1239 Double_t fieldZ = AliTracker::GetBz(); //
1240 AliESDtrack * track = esd->GetTrack(itrack);
1241 if (!track) continue;
1242 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
a28195f7 1243 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
ff8bb5ae 1244 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1245 }
1246
5f8272e1 1247 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1248 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1249
596a855f 1250 return kTRUE;
1251}
1252
1253//_____________________________________________________________________________
24f7a148 1254Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1255{
1256// fill the event summary data
1257
030b532d 1258 TStopwatch stopwatch;
1259 stopwatch.Start();
815c2b38 1260 AliInfo("filling ESD");
030b532d 1261
596a855f 1262 TString detStr = detectors;
b8cd5251 1263 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1264 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1265 AliReconstructor* reconstructor = GetReconstructor(iDet);
1266 if (!reconstructor) continue;
1267
1268 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1269 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 1270 TTree* clustersTree = NULL;
1271 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1272 fLoader[iDet]->LoadRecPoints("read");
1273 clustersTree = fLoader[iDet]->TreeR();
1274 if (!clustersTree) {
1275 AliError(Form("Can't get the %s clusters tree",
1276 fgkDetectorName[iDet]));
1277 if (fStopOnError) return kFALSE;
1278 }
1279 }
1280 if (fRawReader && !reconstructor->HasDigitConversion()) {
1281 reconstructor->FillESD(fRawReader, clustersTree, esd);
1282 } else {
1283 TTree* digitsTree = NULL;
1284 if (fLoader[iDet]) {
1285 fLoader[iDet]->LoadDigits("read");
1286 digitsTree = fLoader[iDet]->TreeD();
1287 if (!digitsTree) {
1288 AliError(Form("Can't get the %s digits tree",
1289 fgkDetectorName[iDet]));
1290 if (fStopOnError) return kFALSE;
1291 }
1292 }
1293 reconstructor->FillESD(digitsTree, clustersTree, esd);
1294 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1295 }
1296 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1297 fLoader[iDet]->UnloadRecPoints();
1298 }
1299
b8cd5251 1300 if (fRawReader) {
1301 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1302 } else {
1303 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1304 }
b8cd5251 1305 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1306 }
1307 }
1308
1309 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1310 AliError(Form("the following detectors were not found: %s",
1311 detStr.Data()));
596a855f 1312 if (fStopOnError) return kFALSE;
1313 }
1314
5f8272e1 1315 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1316 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1317
596a855f 1318 return kTRUE;
1319}
1320
b647652d 1321//_____________________________________________________________________________
1322Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1323{
1324 // Reads the trigger decision which is
1325 // stored in Trigger.root file and fills
1326 // the corresponding esd entries
1327
1328 AliInfo("Filling trigger information into the ESD");
1329
1330 if (fRawReader) {
1331 AliCTPRawStream input(fRawReader);
1332 if (!input.Next()) {
1333 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1334 return kFALSE;
1335 }
1336 esd->SetTriggerMask(input.GetClassMask());
1337 esd->SetTriggerCluster(input.GetClusterMask());
1338 }
1339 else {
1340 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1341 if (runloader) {
1342 if (!runloader->LoadTrigger()) {
1343 AliCentralTrigger *aCTP = runloader->GetTrigger();
1344 esd->SetTriggerMask(aCTP->GetClassMask());
1345 esd->SetTriggerCluster(aCTP->GetClusterMask());
1346 }
1347 else {
1348 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1349 return kFALSE;
1350 }
1351 }
1352 else {
1353 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1354 return kFALSE;
1355 }
1356 }
1357
1358 return kTRUE;
1359}
596a855f 1360
001397cd 1361
1362
1363
1364
1365//_____________________________________________________________________________
1366Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1367{
1368 //
1369 // Filling information from RawReader Header
1370 //
1371
1372 AliInfo("Filling information from RawReader Header");
31fd97b2 1373 esd->SetBunchCrossNumber(0);
1374 esd->SetOrbitNumber(0);
9bcc1e45 1375 esd->SetPeriodNumber(0);
001397cd 1376 esd->SetTimeStamp(0);
1377 esd->SetEventType(0);
1378 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1379 if (eventHeader){
9bcc1e45 1380
1381 const UInt_t *id = eventHeader->GetP("Id");
1382 esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1383 esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1384 esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1385
001397cd 1386 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
31fd97b2 1387 esd->SetEventType((eventHeader->Get("Type")));
001397cd 1388 }
1389
1390 return kTRUE;
1391}
1392
1393
596a855f 1394//_____________________________________________________________________________
1395Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1396{
1397// check whether detName is contained in detectors
1398// if yes, it is removed from detectors
1399
1400 // check if all detectors are selected
1401 if ((detectors.CompareTo("ALL") == 0) ||
1402 detectors.BeginsWith("ALL ") ||
1403 detectors.EndsWith(" ALL") ||
1404 detectors.Contains(" ALL ")) {
1405 detectors = "ALL";
1406 return kTRUE;
1407 }
1408
1409 // search for the given detector
1410 Bool_t result = kFALSE;
1411 if ((detectors.CompareTo(detName) == 0) ||
1412 detectors.BeginsWith(detName+" ") ||
1413 detectors.EndsWith(" "+detName) ||
1414 detectors.Contains(" "+detName+" ")) {
1415 detectors.ReplaceAll(detName, "");
1416 result = kTRUE;
1417 }
1418
1419 // clean up the detectors string
1420 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1421 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1422 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1423
1424 return result;
1425}
e583c30d 1426
f08fc9f5 1427//_____________________________________________________________________________
1428Bool_t AliReconstruction::InitRunLoader()
1429{
1430// get or create the run loader
1431
1432 if (gAlice) delete gAlice;
1433 gAlice = NULL;
1434
b26c3770 1435 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1436 // load all base libraries to get the loader classes
1437 TString libs = gSystem->GetLibraries();
1438 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1439 TString detName = fgkDetectorName[iDet];
1440 if (detName == "HLT") continue;
1441 if (libs.Contains("lib" + detName + "base.so")) continue;
1442 gSystem->Load("lib" + detName + "base.so");
1443 }
f08fc9f5 1444 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1445 if (!fRunLoader) {
1446 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1447 CleanUp();
1448 return kFALSE;
1449 }
b26c3770 1450 fRunLoader->CdGAFile();
1451 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1452 if (fRunLoader->LoadgAlice() == 0) {
1453 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1454 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1455 }
f08fc9f5 1456 }
1457 if (!gAlice && !fRawReader) {
1458 AliError(Form("no gAlice object found in file %s",
1459 fGAliceFileName.Data()));
1460 CleanUp();
1461 return kFALSE;
1462 }
1463
1464 } else { // galice.root does not exist
1465 if (!fRawReader) {
1466 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1467 CleanUp();
1468 return kFALSE;
1469 }
1470 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1471 AliConfig::GetDefaultEventFolderName(),
1472 "recreate");
1473 if (!fRunLoader) {
1474 AliError(Form("could not create run loader in file %s",
1475 fGAliceFileName.Data()));
1476 CleanUp();
1477 return kFALSE;
1478 }
1479 fRunLoader->MakeTree("E");
1480 Int_t iEvent = 0;
1481 while (fRawReader->NextEvent()) {
1482 fRunLoader->SetEventNumber(iEvent);
1483 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1484 iEvent, iEvent);
1485 fRunLoader->MakeTree("H");
1486 fRunLoader->TreeE()->Fill();
1487 iEvent++;
1488 }
1489 fRawReader->RewindEvents();
973388c2 1490 if (fNumberOfEventsPerFile > 0)
1491 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1492 else
1493 fRunLoader->SetNumberOfEventsPerFile(iEvent);
f08fc9f5 1494 fRunLoader->WriteHeader("OVERWRITE");
1495 fRunLoader->CdGAFile();
1496 fRunLoader->Write(0, TObject::kOverwrite);
1497// AliTracker::SetFieldMap(???);
1498 }
1499
1500 return kTRUE;
1501}
1502
c757bafd 1503//_____________________________________________________________________________
b8cd5251 1504AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1505{
f08fc9f5 1506// get the reconstructor object and the loader for a detector
c757bafd 1507
b8cd5251 1508 if (fReconstructor[iDet]) return fReconstructor[iDet];
1509
1510 // load the reconstructor object
1511 TPluginManager* pluginManager = gROOT->GetPluginManager();
1512 TString detName = fgkDetectorName[iDet];
1513 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1514 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1515
1516 if (detName == "HLT") {
1517 if (!gROOT->GetClass("AliLevel3")) {
4aa41877 1518 gSystem->Load("libAliHLTSrc.so");
1519 gSystem->Load("libAliHLTMisc.so");
1520 gSystem->Load("libAliHLTHough.so");
1521 gSystem->Load("libAliHLTComp.so");
b8cd5251 1522 }
1523 }
1524
1525 AliReconstructor* reconstructor = NULL;
1526 // first check if a plugin is defined for the reconstructor
1527 TPluginHandler* pluginHandler =
1528 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1529 // if not, add a plugin for it
1530 if (!pluginHandler) {
b8cd5251 1531 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1532 TString libs = gSystem->GetLibraries();
1533 if (libs.Contains("lib" + detName + "base.so") ||
1534 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1535 pluginManager->AddHandler("AliReconstructor", detName,
1536 recName, detName + "rec", recName + "()");
1537 } else {
1538 pluginManager->AddHandler("AliReconstructor", detName,
1539 recName, detName, recName + "()");
c757bafd 1540 }
b8cd5251 1541 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1542 }
1543 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1544 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1545 }
b8cd5251 1546 if (reconstructor) {
1547 TObject* obj = fOptions.FindObject(detName.Data());
1548 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1549 reconstructor->Init(fRunLoader);
b8cd5251 1550 fReconstructor[iDet] = reconstructor;
1551 }
1552
f08fc9f5 1553 // get or create the loader
1554 if (detName != "HLT") {
1555 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1556 if (!fLoader[iDet]) {
1557 AliConfig::Instance()
1558 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1559 detName, detName);
1560 // first check if a plugin is defined for the loader
bb0901a4 1561 pluginHandler =
f08fc9f5 1562 pluginManager->FindHandler("AliLoader", detName);
1563 // if not, add a plugin for it
1564 if (!pluginHandler) {
1565 TString loaderName = "Ali" + detName + "Loader";
1566 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1567 pluginManager->AddHandler("AliLoader", detName,
1568 loaderName, detName + "base",
1569 loaderName + "(const char*, TFolder*)");
1570 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1571 }
1572 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1573 fLoader[iDet] =
1574 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1575 fRunLoader->GetEventFolder());
1576 }
1577 if (!fLoader[iDet]) { // use default loader
1578 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1579 }
1580 if (!fLoader[iDet]) {
1581 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1582 if (fStopOnError) return NULL;
f08fc9f5 1583 } else {
1584 fRunLoader->AddLoader(fLoader[iDet]);
1585 fRunLoader->CdGAFile();
1586 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1587 fRunLoader->Write(0, TObject::kOverwrite);
1588 }
1589 }
1590 }
1591
b8cd5251 1592 return reconstructor;
c757bafd 1593}
1594
2257f27e 1595//_____________________________________________________________________________
1596Bool_t AliReconstruction::CreateVertexer()
1597{
1598// create the vertexer
1599
b8cd5251 1600 fVertexer = NULL;
1601 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1602 if (itsReconstructor) {
b8cd5251 1603 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1604 }
b8cd5251 1605 if (!fVertexer) {
815c2b38 1606 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1607 if (fStopOnError) return kFALSE;
1608 }
1609
1610 return kTRUE;
1611}
1612
24f7a148 1613//_____________________________________________________________________________
b8cd5251 1614Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1615{
f08fc9f5 1616// create the trackers
24f7a148 1617
b8cd5251 1618 TString detStr = detectors;
1619 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1620 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1621 AliReconstructor* reconstructor = GetReconstructor(iDet);
1622 if (!reconstructor) continue;
1623 TString detName = fgkDetectorName[iDet];
1f46a9ae 1624 if (detName == "HLT") {
1625 fRunHLTTracking = kTRUE;
1626 continue;
1627 }
e66fbafb 1628 if (detName == "MUON") {
1629 fRunMuonTracking = kTRUE;
1630 continue;
1631 }
1632
f08fc9f5 1633
1634 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1635 if (!fTracker[iDet] && (iDet < 7)) {
1636 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1637 if (fStopOnError) return kFALSE;
1638 }
1639 }
1640
24f7a148 1641 return kTRUE;
1642}
1643
e583c30d 1644//_____________________________________________________________________________
b26c3770 1645void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1646{
1647// delete trackers and the run loader and close and delete the file
1648
b8cd5251 1649 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1650 delete fReconstructor[iDet];
1651 fReconstructor[iDet] = NULL;
1652 fLoader[iDet] = NULL;
1653 delete fTracker[iDet];
1654 fTracker[iDet] = NULL;
1655 }
1656 delete fVertexer;
1657 fVertexer = NULL;
9178838a 1658 delete fDiamondProfile;
1659 fDiamondProfile = NULL;
e583c30d 1660
1661 delete fRunLoader;
1662 fRunLoader = NULL;
b649205a 1663 delete fRawReader;
1664 fRawReader = NULL;
e583c30d 1665
1666 if (file) {
1667 file->Close();
1668 delete file;
1669 }
b26c3770 1670
1671 if (fileOld) {
1672 fileOld->Close();
1673 delete fileOld;
1674 gSystem->Unlink("AliESDs.old.root");
1675 }
e583c30d 1676}
24f7a148 1677
1678
1679//_____________________________________________________________________________
1680Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1681{
1682// read the ESD event from a file
1683
1684 if (!esd) return kFALSE;
1685 char fileName[256];
1686 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1687 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1688 if (gSystem->AccessPathName(fileName)) return kFALSE;
1689
f3a97c86 1690 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1691 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1692 TFile* file = TFile::Open(fileName);
1693 if (!file || !file->IsOpen()) {
815c2b38 1694 AliError(Form("opening %s failed", fileName));
24f7a148 1695 delete file;
1696 return kFALSE;
1697 }
1698
1699 gROOT->cd();
1700 delete esd;
1701 esd = (AliESD*) file->Get("ESD");
1702 file->Close();
1703 delete file;
1704 return kTRUE;
1705}
1706
1707//_____________________________________________________________________________
1708void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1709{
1710// write the ESD event to a file
1711
1712 if (!esd) return;
1713 char fileName[256];
1714 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1715 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1716
815c2b38 1717 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1718 TFile* file = TFile::Open(fileName, "recreate");
1719 if (!file || !file->IsOpen()) {
815c2b38 1720 AliError(Form("opening %s failed", fileName));
24f7a148 1721 } else {
1722 esd->Write("ESD");
1723 file->Close();
1724 }
1725 delete file;
1726}
f3a97c86 1727
1728
1729
1730
1731//_____________________________________________________________________________
1732void AliReconstruction::CreateTag(TFile* file)
1733{
1ac7d449 1734 //GRP
1735 Float_t lhcLuminosity = 0.0;
1736 TString lhcState = "test";
1737 UInt_t detectorMask = 0;
1738
2bdb9d38 1739 /////////////
1740 //muon code//
1741 ////////////
56982dd1 1742 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1743 //Variables
56982dd1 1744 Double_t fX,fY,fZ ;
1745 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1746 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1747 Int_t fCharge;
1748 TLorentzVector fEPvector;
1749
1750 Float_t fZVertexCut = 10.0;
1751 Float_t fRhoVertexCut = 2.0;
1752
1753 Float_t fLowPtCut = 1.0;
1754 Float_t fHighPtCut = 3.0;
1755 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1756 ////////////
1757
1758 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1759
089bf903 1760 // Creates the tags for all the events in a given ESD file
4302e20f 1761 Int_t ntrack;
089bf903 1762 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1763 Int_t nPos, nNeg, nNeutr;
1764 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1765 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1766 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1767 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1768 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1769 Int_t fVertexflag;
1387d165 1770 Int_t iRunNumber = 0;
547d75a6 1771 TString fVertexName("default");
4302e20f 1772
f3a97c86 1773 AliRunTag *tag = new AliRunTag();
f3a97c86 1774 AliEventTag *evTag = new AliEventTag();
1775 TTree ttag("T","A Tree with event tags");
2bdb9d38 1776 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1777 btag->SetCompressionLevel(9);
2bdb9d38 1778
f3a97c86 1779 AliInfo(Form("Creating the tags......."));
1780
1781 if (!file || !file->IsOpen()) {
1782 AliError(Form("opening failed"));
1783 delete file;
1784 return ;
2bdb9d38 1785 }
d71aa170 1786 Int_t lastEvent = 0;
46698ae4 1787 TTree *b = (TTree*) file->Get("esdTree");
1788 AliESD *esd = new AliESD();
1789 esd->ReadFromTree(b);
1387d165 1790
d71aa170 1791 b->GetEntry(fFirstEvent);
1387d165 1792 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1793
bb0901a4 1794 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1795 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1796 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 1797 ntrack = 0;
1798 nPos = 0;
1799 nNeg = 0;
1800 nNeutr =0;
1801 nK0s = 0;
1802 nNeutrons = 0;
1803 nPi0s = 0;
1804 nGamas = 0;
1805 nProtons = 0;
1806 nKaons = 0;
1807 nPions = 0;
1808 nMuons = 0;
1809 nElectrons = 0;
1810 nCh1GeV = 0;
1811 nCh3GeV = 0;
1812 nCh10GeV = 0;
1813 nMu1GeV = 0;
1814 nMu3GeV = 0;
1815 nMu10GeV = 0;
1816 nEl1GeV = 0;
1817 nEl3GeV = 0;
1818 nEl10GeV = 0;
1819 maxPt = .0;
1820 meanPt = .0;
1821 totalP = .0;
547d75a6 1822 fVertexflag = 0;
d1a50cb5 1823
2bdb9d38 1824 b->GetEntry(iEventNumber);
1387d165 1825 iRunNumber = esd->GetRunNumber();
1826 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1827 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1828 if (!vertexIn) AliError("ESD has not defined vertex.");
1829 if (vertexIn) fVertexName = vertexIn->GetName();
1830 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1831 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1832 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1833 UInt_t status = esdTrack->GetStatus();
f3a97c86 1834
2bdb9d38 1835 //select only tracks with ITS refit
1836 if ((status&AliESDtrack::kITSrefit)==0) continue;
1837 //select only tracks with TPC refit
1838 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1839
2bdb9d38 1840 //select only tracks with the "combined PID"
1841 if ((status&AliESDtrack::kESDpid)==0) continue;
1842 Double_t p[3];
1843 esdTrack->GetPxPyPz(p);
1844 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1845 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1846 totalP += momentum;
1847 meanPt += fPt;
1848 if(fPt > maxPt) maxPt = fPt;
4302e20f 1849
2bdb9d38 1850 if(esdTrack->GetSign() > 0) {
1851 nPos++;
56982dd1 1852 if(fPt > fLowPtCut) nCh1GeV++;
1853 if(fPt > fHighPtCut) nCh3GeV++;
1854 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1855 }
1856 if(esdTrack->GetSign() < 0) {
1857 nNeg++;
56982dd1 1858 if(fPt > fLowPtCut) nCh1GeV++;
1859 if(fPt > fHighPtCut) nCh3GeV++;
1860 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1861 }
1862 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1863
2bdb9d38 1864 //PID
1865 Double_t prob[5];
1866 esdTrack->GetESDpid(prob);
4302e20f 1867
2bdb9d38 1868 Double_t rcc = 0.0;
1869 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1870 if(rcc == 0.0) continue;
1871 //Bayes' formula
1872 Double_t w[5];
1873 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1874
2bdb9d38 1875 //protons
1876 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1877 //kaons
1878 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1879 //pions
1880 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1881 //electrons
1882 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1883 nElectrons++;
56982dd1 1884 if(fPt > fLowPtCut) nEl1GeV++;
1885 if(fPt > fHighPtCut) nEl3GeV++;
1886 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1887 }
1888 ntrack++;
1889 }//track loop
1890
1891 /////////////
1892 //muon code//
1893 ////////////
1894 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1895 // loop over all reconstructed tracks (also first track of combination)
1896 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1897 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1898 if (muonTrack == 0x0) continue;
4302e20f 1899
2bdb9d38 1900 // Coordinates at vertex
56982dd1 1901 fZ = muonTrack->GetZ();
1902 fY = muonTrack->GetBendingCoor();
1903 fX = muonTrack->GetNonBendingCoor();
4302e20f 1904
56982dd1 1905 fThetaX = muonTrack->GetThetaX();
1906 fThetaY = muonTrack->GetThetaY();
4302e20f 1907
56982dd1 1908 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1909 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1910 fPxRec = fPzRec * TMath::Tan(fThetaX);
1911 fPyRec = fPzRec * TMath::Tan(fThetaY);
1912 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1913
2bdb9d38 1914 //ChiSquare of the track if needed
56982dd1 1915 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1916 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1917 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1918
2bdb9d38 1919 // total number of muons inside a vertex cut
56982dd1 1920 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1921 nMuons++;
56982dd1 1922 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1923 nMu1GeV++;
56982dd1 1924 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1925 nMu3GeV++;
56982dd1 1926 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1927 nMu10GeV++;
1928 }
1929 }
1930 }
1931 }
1932 }//muon track loop
1933
1934 // Fill the event tags
1935 if(ntrack != 0)
1936 meanPt = meanPt/ntrack;
1937
1938 evTag->SetEventId(iEventNumber+1);
547d75a6 1939 if (vertexIn) {
1940 evTag->SetVertexX(vertexIn->GetXv());
1941 evTag->SetVertexY(vertexIn->GetYv());
1942 evTag->SetVertexZ(vertexIn->GetZv());
1943 evTag->SetVertexZError(vertexIn->GetZRes());
1944 }
d1a50cb5 1945 evTag->SetVertexFlag(fVertexflag);
1946
2bdb9d38 1947 evTag->SetT0VertexZ(esd->GetT0zVertex());
1948
8bd8ac26 1949 evTag->SetTriggerMask(esd->GetTriggerMask());
1950 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1951
32a5cab4 1952 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1953 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1954 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1955 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1956 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1957 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1958
1959
1960 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1961 evTag->SetNumOfPosTracks(nPos);
1962 evTag->SetNumOfNegTracks(nNeg);
1963 evTag->SetNumOfNeutrTracks(nNeutr);
1964
1965 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1966 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1967 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1968 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1969
1970 evTag->SetNumOfProtons(nProtons);
1971 evTag->SetNumOfKaons(nKaons);
1972 evTag->SetNumOfPions(nPions);
1973 evTag->SetNumOfMuons(nMuons);
1974 evTag->SetNumOfElectrons(nElectrons);
1975 evTag->SetNumOfPhotons(nGamas);
1976 evTag->SetNumOfPi0s(nPi0s);
1977 evTag->SetNumOfNeutrons(nNeutrons);
1978 evTag->SetNumOfKaon0s(nK0s);
1979
1980 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1981 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1982 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1983 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1984 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1985 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1986 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1987 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1988 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1989
85c60a8e 1990 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1991 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1992
1993 evTag->SetTotalMomentum(totalP);
1994 evTag->SetMeanPt(meanPt);
1995 evTag->SetMaxPt(maxPt);
1996
1ac7d449 1997 tag->SetLHCTag(lhcLuminosity,lhcState);
1998 tag->SetDetectorTag(detectorMask);
1999
1387d165 2000 tag->SetRunId(iInitRunNumber);
2bdb9d38 2001 tag->AddEventTag(*evTag);
2002 }
bb0901a4 2003 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 2004 else lastEvent = fLastEvent;
f3a97c86 2005
2006 ttag.Fill();
2007 tag->Clear();
2008
2009 char fileName[256];
2010 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 2011 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 2012 AliInfo(Form("writing tags to file %s", fileName));
2013 AliDebug(1, Form("writing tags to file %s", fileName));
2014
2015 TFile* ftag = TFile::Open(fileName, "recreate");
2016 ftag->cd();
2017 ttag.Write();
2018 ftag->Close();
2019 file->cd();
2020 delete tag;
f3a97c86 2021 delete evTag;
2022}
2023
a7807689 2024//_____________________________________________________________________________
f29f1726 2025void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
a7807689 2026{
f29f1726 2027 // write all files from the given esd file to an aod file
2028
2029 // create an AliAOD object
2030 AliAODEvent *aod = new AliAODEvent();
2031 aod->CreateStdContent();
2032
2033 // go to the file
2034 aodFile->cd();
2035
2036 // create the tree
2037 TTree *aodTree = new TTree("AOD", "AliAOD tree");
2038 aodTree->Branch(aod->GetList());
2039
2040 // connect to ESD
2041 TTree *t = (TTree*) esdFile->Get("esdTree");
2042 TBranch *b = t->GetBranch("ESD");
2043 AliESD *esd = 0;
2044 b->SetAddress(&esd);
2045
2046 Int_t nEvents = b->GetEntries();
2047
2048 // set arrays and pointers
2049 Float_t posF[3];
2050 Double_t pos[3];
2051 Double_t p[3];
2052 Double_t covVtx[6];
2053 Double_t covTr[21];
2054 Double_t pid[10];
2055
2056 // loop over events and fill them
2057 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2058 b->GetEntry(iEvent);
2059
2060 // Multiplicity information needed by the header (to be revised!)
2061 Int_t nTracks = esd->GetNumberOfTracks();
2062 Int_t nPosTracks = 0;
2063 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2064 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2065
2066 // create the header
2067 aod->AddHeader(new AliAODHeader(esd->GetRunNumber(),
2068 esd->GetBunchCrossNumber(),
2069 esd->GetOrbitNumber(),
2070 esd->GetPeriodNumber(),
2071 nTracks,
2072 nPosTracks,
2073 nTracks-nPosTracks,
2074 esd->GetMagneticField(),
2075 -999., // fill muon magnetic field
2076 -999., // centrality; to be filled, still
2077 esd->GetZDCN1Energy(),
2078 esd->GetZDCP1Energy(),
2079 esd->GetZDCN2Energy(),
2080 esd->GetZDCP2Energy(),
2081 esd->GetZDCEMEnergy(),
2082 esd->GetTriggerMask(),
2083 esd->GetTriggerCluster(),
2084 esd->GetEventType()));
2085
2086 Int_t nV0s = esd->GetNumberOfV0s();
2087 Int_t nCascades = esd->GetNumberOfCascades();
2088 Int_t nKinks = esd->GetNumberOfKinks();
2089 Int_t nVertices = nV0s + nCascades + nKinks;
2090
2091 aod->ResetStd(nTracks, nVertices);
2092 AliAODTrack *aodTrack;
2093
2094
2095 // Array to take into account the tracks already added to the AOD
2096 Bool_t * usedTrack = NULL;
2097 if (nTracks>0) {
2098 usedTrack = new Bool_t[nTracks];
2099 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2100 }
2101 // Array to take into account the V0s already added to the AOD
2102 Bool_t * usedV0 = NULL;
2103 if (nV0s>0) {
2104 usedV0 = new Bool_t[nV0s];
2105 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2106 }
2107 // Array to take into account the kinks already added to the AOD
2108 Bool_t * usedKink = NULL;
2109 if (nKinks>0) {
2110 usedKink = new Bool_t[nKinks];
2111 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2112 }
2113
2114 // Access to the AOD container of vertices
2115 TClonesArray &vertices = *(aod->GetVertices());
2116 Int_t jVertices=0;
2117
2118 // Access to the AOD container of tracks
2119 TClonesArray &tracks = *(aod->GetTracks());
2120 Int_t jTracks=0;
2121
2122 // Add primary vertex. The primary tracks will be defined
2123 // after the loops on the composite objects (V0, cascades, kinks)
2124 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2125
2126 vtx->GetXYZ(pos); // position
2127 vtx->GetCovMatrix(covVtx); //covariance matrix
2128
2129 AliAODVertex * primary = new(vertices[jVertices++])
2130 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2131
2132 // Create vertices starting from the most complex objects
2133
2134 // Cascades
2135 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2136 AliESDcascade *cascade = esd->GetCascade(nCascade);
2137
2138 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2139 cascade->GetPosCovXi(covVtx);
2140
2141 // Add the cascade vertex
2142 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2143 covVtx,
2144 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2145 primary,
2146 AliAODVertex::kCascade);
2147
2148 primary->AddDaughter(vcascade);
2149
2150 // Add the V0 from the cascade. The ESD class have to be optimized...
2151 // Now we have to search for the corresponding Vo in the list of V0s
2152 // using the indeces of the positive and negative tracks
2153
2154 Int_t posFromV0 = cascade->GetPindex();
2155 Int_t negFromV0 = cascade->GetNindex();
2156
2157
2158 AliESDv0 * v0 = 0x0;
2159 Int_t indV0 = -1;
2160
2161 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2162
2163 v0 = esd->GetV0(iV0);
2164 Int_t posV0 = v0->GetPindex();
2165 Int_t negV0 = v0->GetNindex();
2166
2167 if (posV0==posFromV0 && negV0==negFromV0) {
2168 indV0 = iV0;
2169 break;
2170 }
2171 }
2172
2173 AliAODVertex * vV0FromCascade = 0x0;
2174
2175 if (indV0>-1 && !usedV0[indV0] ) {
2176
2177 // the V0 exists in the array of V0s and is not used
2178
2179 usedV0[indV0] = kTRUE;
2180
2181 v0->GetXYZ(pos[0], pos[1], pos[2]);
2182 v0->GetPosCov(covVtx);
2183
2184 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2185 covVtx,
2186 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2187 vcascade,
2188 AliAODVertex::kV0);
2189 } else {
2190
2191 // the V0 doesn't exist in the array of V0s or was used
2192 cerr << "Error: event " << iEvent << " cascade " << nCascade
2193 << " The V0 " << indV0
2194 << " doesn't exist in the array of V0s or was used!" << endl;
2195
2196 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2197 cascade->GetPosCov(covVtx);
2198
2199 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2200 covVtx,
2201 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2202 vcascade,
2203 AliAODVertex::kV0);
2204 vcascade->AddDaughter(vV0FromCascade);
2205 }
2206
2207 // Add the positive tracks from the V0
2208
2209 if (! usedTrack[posFromV0]) {
2210
2211 usedTrack[posFromV0] = kTRUE;
2212
2213 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2214 esdTrack->GetPxPyPz(p);
2215 esdTrack->GetXYZ(pos);
2216 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2217 esdTrack->GetESDpid(pid);
2218
2219 vV0FromCascade->AddDaughter(aodTrack =
2220 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2221 esdTrack->GetLabel(),
2222 p,
2223 kTRUE,
2224 pos,
2225 kFALSE,
2226 covTr,
2227 (Short_t)esdTrack->GetSign(),
2228 esdTrack->GetITSClusterMap(),
2229 pid,
2230 vV0FromCascade,
2231 kTRUE, // check if this is right
2232 kFALSE, // check if this is right
2233 AliAODTrack::kSecondary)
2234 );
2235 aodTrack->ConvertAliPIDtoAODPID();
2236 }
2237 else {
2238 cerr << "Error: event " << iEvent << " cascade " << nCascade
2239 << " track " << posFromV0 << " has already been used!" << endl;
2240 }
2241
2242 // Add the negative tracks from the V0
2243
2244 if (!usedTrack[negFromV0]) {
2245
2246 usedTrack[negFromV0] = kTRUE;
2247
2248 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2249 esdTrack->GetPxPyPz(p);
2250 esdTrack->GetXYZ(pos);
2251 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2252 esdTrack->GetESDpid(pid);
2253
2254 vV0FromCascade->AddDaughter(aodTrack =
2255 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2256 esdTrack->GetLabel(),
2257 p,
2258 kTRUE,
2259 pos,
2260 kFALSE,
2261 covTr,
2262 (Short_t)esdTrack->GetSign(),
2263 esdTrack->GetITSClusterMap(),
2264 pid,
2265 vV0FromCascade,
2266 kTRUE, // check if this is right
2267 kFALSE, // check if this is right
2268 AliAODTrack::kSecondary)
2269 );
2270 aodTrack->ConvertAliPIDtoAODPID();
2271 }
2272 else {
2273 cerr << "Error: event " << iEvent << " cascade " << nCascade
2274 << " track " << negFromV0 << " has already been used!" << endl;
2275 }
2276
2277 // Add the bachelor track from the cascade
2278
2279 Int_t bachelor = cascade->GetBindex();
2280
2281 if(!usedTrack[bachelor]) {
2282
2283 usedTrack[bachelor] = kTRUE;
2284
2285 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2286 esdTrack->GetPxPyPz(p);
2287 esdTrack->GetXYZ(pos);
2288 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2289 esdTrack->GetESDpid(pid);
2290
2291 vcascade->AddDaughter(aodTrack =
2292 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2293 esdTrack->GetLabel(),
2294 p,
2295 kTRUE,
2296 pos,
2297 kFALSE,
2298 covTr,
2299 (Short_t)esdTrack->GetSign(),
2300 esdTrack->GetITSClusterMap(),
2301 pid,
2302 vcascade,
2303 kTRUE, // check if this is right
2304 kFALSE, // check if this is right
2305 AliAODTrack::kSecondary)
2306 );
2307 aodTrack->ConvertAliPIDtoAODPID();
2308 }
2309 else {
2310 cerr << "Error: event " << iEvent << " cascade " << nCascade
2311 << " track " << bachelor << " has already been used!" << endl;
2312 }
2313
2314 // Add the primary track of the cascade (if any)
2315
2316 } // end of the loop on cascades
2317
2318 // V0s
2319
2320 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2321
2322 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2323
2324 AliESDv0 *v0 = esd->GetV0(nV0);
2325
2326 v0->GetXYZ(pos[0], pos[1], pos[2]);
2327 v0->GetPosCov(covVtx);
2328
2329 AliAODVertex * vV0 =
2330 new(vertices[jVertices++]) AliAODVertex(pos,
2331 covVtx,
2332 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2333 primary,
2334 AliAODVertex::kV0);
2335 primary->AddDaughter(vV0);
2336
2337 Int_t posFromV0 = v0->GetPindex();
2338 Int_t negFromV0 = v0->GetNindex();
2339
2340 // Add the positive tracks from the V0
2341
2342 if (!usedTrack[posFromV0]) {
2343
2344 usedTrack[posFromV0] = kTRUE;
2345
2346 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2347 esdTrack->GetPxPyPz(p);
2348 esdTrack->GetXYZ(pos);
2349 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2350 esdTrack->GetESDpid(pid);
2351
2352 vV0->AddDaughter(aodTrack =
2353 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2354 esdTrack->GetLabel(),
2355 p,
2356 kTRUE,
2357 pos,
2358 kFALSE,
2359 covTr,
2360 (Short_t)esdTrack->GetSign(),
2361 esdTrack->GetITSClusterMap(),
2362 pid,
2363 vV0,
2364 kTRUE, // check if this is right
2365 kFALSE, // check if this is right
2366 AliAODTrack::kSecondary)
2367 );
2368 aodTrack->ConvertAliPIDtoAODPID();
2369 }
2370 else {
2371 cerr << "Error: event " << iEvent << " V0 " << nV0
2372 << " track " << posFromV0 << " has already been used!" << endl;
2373 }
a7807689 2374
f29f1726 2375 // Add the negative tracks from the V0
2376
2377 if (!usedTrack[negFromV0]) {
2378
2379 usedTrack[negFromV0] = kTRUE;
2380
2381 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2382 esdTrack->GetPxPyPz(p);
2383 esdTrack->GetXYZ(pos);
2384 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2385 esdTrack->GetESDpid(pid);
2386
2387 vV0->AddDaughter(aodTrack =
2388 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2389 esdTrack->GetLabel(),
2390 p,
2391 kTRUE,
2392 pos,
2393 kFALSE,
2394 covTr,
2395 (Short_t)esdTrack->GetSign(),
2396 esdTrack->GetITSClusterMap(),
2397 pid,
2398 vV0,
2399 kTRUE, // check if this is right
2400 kFALSE, // check if this is right
2401 AliAODTrack::kSecondary)
2402 );
2403 aodTrack->ConvertAliPIDtoAODPID();
2404 }
2405 else {
2406 cerr << "Error: event " << iEvent << " V0 " << nV0
2407 << " track " << negFromV0 << " has already been used!" << endl;
2408 }
2409
2410 } // end of the loop on V0s
2411
2412 // Kinks: it is a big mess the access to the information in the kinks
2413 // The loop is on the tracks in order to find the mother and daugther of each kink
2414
2415
2416 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2417
2418
2419 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2420
2421 Int_t ikink = esdTrack->GetKinkIndex(0);
2422
2423 if (ikink) {
2424 // Negative kink index: mother, positive: daughter
2425
2426 // Search for the second track of the kink
2427
2428 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2429
2430 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2431
2432 Int_t jkink = esdTrack1->GetKinkIndex(0);
2433
2434 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2435
2436 // The two tracks are from the same kink
2437
2438 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2439
2440 Int_t imother = -1;
2441 Int_t idaughter = -1;
2442
2443 if (ikink<0 && jkink>0) {
2444
2445 imother = iTrack;
2446 idaughter = jTrack;
2447 }
2448 else if (ikink>0 && jkink<0) {
2449
2450 imother = jTrack;
2451 idaughter = iTrack;
2452 }
2453 else {
2454 cerr << "Error: Wrong combination of kink indexes: "
2455 << ikink << " " << jkink << endl;
2456 continue;
2457 }
2458
2459 // Add the mother track
2460
2461 AliAODTrack * mother = NULL;
2462
2463 if (!usedTrack[imother]) {
2464
2465 usedTrack[imother] = kTRUE;
2466
2467 AliESDtrack *esdTrack = esd->GetTrack(imother);
2468 esdTrack->GetPxPyPz(p);
2469 esdTrack->GetXYZ(pos);
2470 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2471 esdTrack->GetESDpid(pid);
2472
2473 mother =
2474 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2475 esdTrack->GetLabel(),
2476 p,
2477 kTRUE,
2478 pos,
2479 kFALSE,
2480 covTr,
2481 (Short_t)esdTrack->GetSign(),
2482 esdTrack->GetITSClusterMap(),
2483 pid,
2484 primary,
2485 kTRUE, // check if this is right
2486 kTRUE, // check if this is right
2487 AliAODTrack::kPrimary);
2488 primary->AddDaughter(mother);
2489 mother->ConvertAliPIDtoAODPID();
2490 }
2491 else {
2492 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2493 << " track " << imother << " has already been used!" << endl;
2494 }
2495
2496 // Add the kink vertex
2497 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2498
2499 AliAODVertex * vkink =
2500 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2501 NULL,
2502 0.,
2503 mother,
2504 AliAODVertex::kKink);
2505 // Add the daughter track
2506
2507 AliAODTrack * daughter = NULL;
2508
2509 if (!usedTrack[idaughter]) {
2510
2511 usedTrack[idaughter] = kTRUE;
2512
2513 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2514 esdTrack->GetPxPyPz(p);
2515 esdTrack->GetXYZ(pos);
2516 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2517 esdTrack->GetESDpid(pid);
2518
2519 daughter =
2520 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2521 esdTrack->GetLabel(),
2522 p,
2523 kTRUE,
2524 pos,
2525 kFALSE,
2526 covTr,
2527 (Short_t)esdTrack->GetSign(),
2528 esdTrack->GetITSClusterMap(),
2529 pid,
2530 vkink,
2531 kTRUE, // check if this is right
2532 kTRUE, // check if this is right
2533 AliAODTrack::kPrimary);
2534 vkink->AddDaughter(daughter);
2535 daughter->ConvertAliPIDtoAODPID();
2536 }
2537 else {
2538 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2539 << " track " << idaughter << " has already been used!" << endl;
2540 }
2541
2542
2543 }
2544 }
2545
2546 }
2547
2548 }
2549
2550
2551 // Tracks (primary and orphan)
2552
2553 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2554
2555
2556 if (usedTrack[nTrack]) continue;
2557
2558 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2559 esdTrack->GetPxPyPz(p);
2560 esdTrack->GetXYZ(pos);
2561 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2562 esdTrack->GetESDpid(pid);
2563
2564 Float_t impactXY, impactZ;
2565
2566 esdTrack->GetImpactParameters(impactXY,impactZ);
2567
2568 if (impactXY<3) {
2569 // track inside the beam pipe
2570
2571 primary->AddDaughter(aodTrack =
2572 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2573 esdTrack->GetLabel(),
2574 p,
2575 kTRUE,
2576 pos,
2577 kFALSE,
2578 covTr,
2579 (Short_t)esdTrack->GetSign(),
2580 esdTrack->GetITSClusterMap(),
2581 pid,
2582 primary,
2583 kTRUE, // check if this is right
2584 kTRUE, // check if this is right
2585 AliAODTrack::kPrimary)
2586 );
2587 aodTrack->ConvertAliPIDtoAODPID();
2588 }
2589 else {
2590 // outside the beam pipe: orphan track
2591 aodTrack =
2592 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2593 esdTrack->GetLabel(),
2594 p,
2595 kTRUE,
2596 pos,
2597 kFALSE,
2598 covTr,
2599 (Short_t)esdTrack->GetSign(),
2600 esdTrack->GetITSClusterMap(),
2601 pid,
2602 NULL,
2603 kFALSE, // check if this is right
2604 kFALSE, // check if this is right
2605 AliAODTrack::kOrphan);
2606 aodTrack->ConvertAliPIDtoAODPID();
2607 }
2608 } // end of loop on tracks
2609
2610 // muon tracks
2611 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2612 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2613
2614 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2615 p[0] = esdMuTrack->Px();
2616 p[1] = esdMuTrack->Py();
2617 p[2] = esdMuTrack->Pz();
2618 pos[0] = primary->GetX();
2619 pos[1] = primary->GetY();
2620 pos[2] = primary->GetZ();
2621
2622 // has to be changed once the muon pid is provided by the ESD
2623 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2624
2625 primary->AddDaughter(
2626 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2627 0, // no label provided
2628 p,
2629 kTRUE,
2630 pos,
2631 kFALSE,
2632 NULL, // no covariance matrix provided
2633 (Short_t)-99, // no charge provided
2634 0, // no ITSClusterMap
2635 pid,
2636 primary,
2637 kTRUE, // check if this is right
2638 kTRUE, // not used for vertex fit
2639 AliAODTrack::kPrimary)
2640 );
2641 }
2642
2643 // Access to the AOD container of clusters
2644 TClonesArray &clusters = *(aod->GetClusters());
2645 Int_t jClusters=0;
2646
2647 // Calo Clusters
2648 Int_t nClusters = esd->GetNumberOfCaloClusters();
2649
2650 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2651
2652 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2653
2654 Int_t id = cluster->GetID();
2655 Int_t label = -1;
2656 Float_t energy = cluster->GetClusterEnergy();
2657 cluster->GetGlobalPosition(posF);
2658 AliAODVertex *prodVertex = primary;
2659 AliAODTrack *primTrack = NULL;
2660 Char_t ttype=AliAODCluster::kUndef;
2661
2662 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2663 else if (cluster->IsEMCAL()) {
2664
2665 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2666 ttype = AliAODCluster::kEMCALPseudoCluster;
2667 else
2668 ttype = AliAODCluster::kEMCALClusterv1;
2669
2670 }
2671
2672 new(clusters[jClusters++]) AliAODCluster(id,
2673 label,
2674 energy,
2675 pos,
2676 NULL, // no covariance matrix provided
2677 NULL, // no pid for clusters provided
2678 prodVertex,
2679 primTrack,
2680 ttype);
2681
2682 } // end of loop on calo clusters
2683
2684 delete [] usedTrack;
2685 delete [] usedV0;
2686 delete [] usedKink;
2687
2688 // fill the tree for this event
2689 aodTree->Fill();
2690 } // end of event loop
2691
2692 aodTree->GetUserInfo()->Add(aod);
2693
2694 // close ESD file
2695 esdFile->Close();
2696
2697 // write the tree to the specified file
2698 aodFile = aodTree->GetCurrentFile();
2699 aodFile->cd();
2700 aodTree->Write();
2701
a7807689 2702 return;
2703}
2704
2705
98937d93 2706void AliReconstruction::WriteAlignmentData(AliESD* esd)
2707{
2708 // Write space-points which are then used in the alignment procedures
2709 // For the moment only ITS, TRD and TPC
2710
2711 // Load TOF clusters
d528ee75 2712 if (fTracker[3]){
2713 fLoader[3]->LoadRecPoints("read");
2714 TTree* tree = fLoader[3]->TreeR();
2715 if (!tree) {
2716 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2717 return;
2718 }
2719 fTracker[3]->LoadClusters(tree);
98937d93 2720 }
98937d93 2721 Int_t ntracks = esd->GetNumberOfTracks();
2722 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2723 {
2724 AliESDtrack *track = esd->GetTrack(itrack);
2725 Int_t nsp = 0;
ef7253ac 2726 Int_t idx[200];
98937d93 2727 for (Int_t iDet = 3; iDet >= 0; iDet--)
2728 nsp += track->GetNcls(iDet);
2729 if (nsp) {
2730 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2731 track->SetTrackPointArray(sp);
2732 Int_t isptrack = 0;
2733 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2734 AliTracker *tracker = fTracker[iDet];
2735 if (!tracker) continue;
2736 Int_t nspdet = track->GetNcls(iDet);
98937d93 2737 if (nspdet <= 0) continue;
2738 track->GetClusters(iDet,idx);
2739 AliTrackPoint p;
2740 Int_t isp = 0;
2741 Int_t isp2 = 0;
2742 while (isp < nspdet) {
2743 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 2744 const Int_t kNTPCmax = 159;
2745 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2746 if (!isvalid) continue;
2747 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2748 }
98937d93 2749 }
2750 }
2751 }
d528ee75 2752 if (fTracker[3]){
2753 fTracker[3]->UnloadClusters();
2754 fLoader[3]->UnloadRecPoints();
2755 }
98937d93 2756}
2e3550da 2757
2758//_____________________________________________________________________________
2759void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2760{
2761 // The method reads the raw-data error log
2762 // accumulated within the rawReader.
2763 // It extracts the raw-data errors related to
2764 // the current event and stores them into
2765 // a TClonesArray inside the esd object.
2766
2767 if (!fRawReader) return;
2768
2769 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2770
2771 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2772 if (!log) continue;
2773 if (iEvent != log->GetEventNumber()) continue;
2774
2775 esd->AddRawDataErrorLog(log);
2776 }
2777
2778}
46698ae4 2779
2780TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2781 ifstream in;
2782 in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2783 Int_t kBytes = (Int_t)in.tellg();
2784 printf("Size: %d \n",kBytes);
2785 TNamed *fn = 0;
2786 if(in.good()){
2787 char* memblock = new char [kBytes];
2788 in.seekg (0, ios::beg);
2789 in.read (memblock, kBytes);
2790 in.close();
2791 TString fData(memblock,kBytes);
2792 fn = new TNamed(fName,fData);
2793 printf("fData Size: %d \n",fData.Sizeof());
2794 printf("fName Size: %d \n",fName.Sizeof());
2795 printf("fn Size: %d \n",fn->Sizeof());
2796 delete[] memblock;
2797 }
2798 else{
2799 AliInfo(Form("Could not Open %s\n",fPath.Data()));
2800 }
2801
2802 return fn;
2803}
2804
2805void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2806
2807 // This is not really needed in AliReconstruction at the moment
2808 // but can serve as a template
2809
2810 TList *fList = fTree->GetUserInfo();
2811 TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2812 printf("fn Size: %d \n",fn->Sizeof());
2813
2814 TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2815 const char* cdata = fn->GetTitle();
2816 printf("fTmp Size %d\n",fTmp.Sizeof());
2817
2818 int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2819 printf("calculated size %d\n",size);
2820 ofstream out(fName.Data(),ios::out | ios::binary);
2821 out.write(cdata,size);
2822 out.close();
2823
2824}
2825