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