Adding Print method (Laurent)
[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
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
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
426//_____________________________________________________________________________
6bae477a 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
883//_____________________________________________________________________________
b26c3770 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
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
1017//_____________________________________________________________________________
1f46a9ae 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
1077//_____________________________________________________________________________
e66fbafb 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
1133//_____________________________________________________________________________
2257f27e 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
1435//_____________________________________________________________________________
f08fc9f5 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
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
1603//_____________________________________________________________________________
2257f27e 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
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
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");
53ec9628 2050 AliESD *esd = new AliESD();
2051 esd->ReadFromTree(t);
f29f1726 2052
53ec9628 2053 Int_t nEvents = t->GetEntries();
f29f1726 2054
2055 // set arrays and pointers
2056 Float_t posF[3];
2057 Double_t pos[3];
2058 Double_t p[3];
2059 Double_t covVtx[6];
2060 Double_t covTr[21];
2061 Double_t pid[10];
2062
2063 // loop over events and fill them
2064 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
53ec9628 2065 t->GetEntry(iEvent);
f29f1726 2066
2067 // Multiplicity information needed by the header (to be revised!)
2068 Int_t nTracks = esd->GetNumberOfTracks();
2069 Int_t nPosTracks = 0;
2070 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2071 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2072
a1d4139d 2073 // Update the header
2074 AliAODHeader* header = aod->GetHeader();
2075
ade23daf 2076 header->SetRunNumber (esd->GetRunNumber() );
2077 header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
2078 header->SetOrbitNumber (esd->GetOrbitNumber() );
2079 header->SetPeriodNumber (esd->GetPeriodNumber() );
2080 header->SetTriggerMask (esd->GetTriggerMask() );
2081 header->SetTriggerCluster (esd->GetTriggerCluster() );
2082 header->SetEventType (esd->GetEventType() );
2083 header->SetMagneticField (esd->GetMagneticField() );
2084 header->SetZDCN1Energy (esd->GetZDCN1Energy() );
2085 header->SetZDCP1Energy (esd->GetZDCP1Energy() );
2086 header->SetZDCN2Energy (esd->GetZDCN2Energy() );
2087 header->SetZDCP2Energy (esd->GetZDCP2Energy() );
2088 header->SetZDCEMEnergy (esd->GetZDCEMEnergy() );
a1d4139d 2089 header->SetRefMultiplicity (nTracks);
2090 header->SetRefMultiplicityPos(nPosTracks);
2091 header->SetRefMultiplicityNeg(nTracks - nPosTracks);
2092 header->SetMuonMagFieldScale(-999.); // FIXME
2093 header->SetCentrality(-999.); // FIXME
2094//
2095//
f29f1726 2096
2097 Int_t nV0s = esd->GetNumberOfV0s();
2098 Int_t nCascades = esd->GetNumberOfCascades();
2099 Int_t nKinks = esd->GetNumberOfKinks();
2100 Int_t nVertices = nV0s + nCascades + nKinks;
2101
2102 aod->ResetStd(nTracks, nVertices);
2103 AliAODTrack *aodTrack;
2104
2105
2106 // Array to take into account the tracks already added to the AOD
2107 Bool_t * usedTrack = NULL;
2108 if (nTracks>0) {
2109 usedTrack = new Bool_t[nTracks];
2110 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2111 }
2112 // Array to take into account the V0s already added to the AOD
2113 Bool_t * usedV0 = NULL;
2114 if (nV0s>0) {
2115 usedV0 = new Bool_t[nV0s];
2116 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2117 }
2118 // Array to take into account the kinks already added to the AOD
2119 Bool_t * usedKink = NULL;
2120 if (nKinks>0) {
2121 usedKink = new Bool_t[nKinks];
2122 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2123 }
2124
2125 // Access to the AOD container of vertices
2126 TClonesArray &vertices = *(aod->GetVertices());
2127 Int_t jVertices=0;
2128
2129 // Access to the AOD container of tracks
2130 TClonesArray &tracks = *(aod->GetTracks());
2131 Int_t jTracks=0;
2132
2133 // Add primary vertex. The primary tracks will be defined
2134 // after the loops on the composite objects (V0, cascades, kinks)
2135 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2136
2137 vtx->GetXYZ(pos); // position
2138 vtx->GetCovMatrix(covVtx); //covariance matrix
2139
2140 AliAODVertex * primary = new(vertices[jVertices++])
2141 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2142
2143 // Create vertices starting from the most complex objects
2144
2145 // Cascades
2146 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2147 AliESDcascade *cascade = esd->GetCascade(nCascade);
2148
2149 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2150 cascade->GetPosCovXi(covVtx);
2151
2152 // Add the cascade vertex
2153 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2154 covVtx,
2155 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2156 primary,
2157 AliAODVertex::kCascade);
2158
2159 primary->AddDaughter(vcascade);
2160
2161 // Add the V0 from the cascade. The ESD class have to be optimized...
2162 // Now we have to search for the corresponding Vo in the list of V0s
2163 // using the indeces of the positive and negative tracks
2164
2165 Int_t posFromV0 = cascade->GetPindex();
2166 Int_t negFromV0 = cascade->GetNindex();
2167
2168
2169 AliESDv0 * v0 = 0x0;
2170 Int_t indV0 = -1;
2171
2172 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2173
2174 v0 = esd->GetV0(iV0);
2175 Int_t posV0 = v0->GetPindex();
2176 Int_t negV0 = v0->GetNindex();
2177
2178 if (posV0==posFromV0 && negV0==negFromV0) {
2179 indV0 = iV0;
2180 break;
2181 }
2182 }
2183
2184 AliAODVertex * vV0FromCascade = 0x0;
2185
2186 if (indV0>-1 && !usedV0[indV0] ) {
2187
2188 // the V0 exists in the array of V0s and is not used
2189
2190 usedV0[indV0] = kTRUE;
2191
2192 v0->GetXYZ(pos[0], pos[1], pos[2]);
2193 v0->GetPosCov(covVtx);
2194
2195 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2196 covVtx,
2197 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2198 vcascade,
2199 AliAODVertex::kV0);
2200 } else {
2201
2202 // the V0 doesn't exist in the array of V0s or was used
2203 cerr << "Error: event " << iEvent << " cascade " << nCascade
2204 << " The V0 " << indV0
2205 << " doesn't exist in the array of V0s or was used!" << endl;
2206
2207 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2208 cascade->GetPosCov(covVtx);
2209
2210 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2211 covVtx,
2212 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2213 vcascade,
2214 AliAODVertex::kV0);
2215 vcascade->AddDaughter(vV0FromCascade);
2216 }
2217
2218 // Add the positive tracks from the V0
2219
2220 if (! usedTrack[posFromV0]) {
2221
2222 usedTrack[posFromV0] = kTRUE;
2223
2224 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2225 esdTrack->GetPxPyPz(p);
2226 esdTrack->GetXYZ(pos);
2227 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2228 esdTrack->GetESDpid(pid);
2229
2230 vV0FromCascade->AddDaughter(aodTrack =
2231 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2232 esdTrack->GetLabel(),
2233 p,
2234 kTRUE,
2235 pos,
2236 kFALSE,
2237 covTr,
2238 (Short_t)esdTrack->GetSign(),
2239 esdTrack->GetITSClusterMap(),
2240 pid,
2241 vV0FromCascade,
2242 kTRUE, // check if this is right
2243 kFALSE, // check if this is right
2244 AliAODTrack::kSecondary)
2245 );
2246 aodTrack->ConvertAliPIDtoAODPID();
2247 }
2248 else {
2249 cerr << "Error: event " << iEvent << " cascade " << nCascade
2250 << " track " << posFromV0 << " has already been used!" << endl;
2251 }
2252
2253 // Add the negative tracks from the V0
2254
2255 if (!usedTrack[negFromV0]) {
2256
2257 usedTrack[negFromV0] = kTRUE;
2258
2259 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2260 esdTrack->GetPxPyPz(p);
2261 esdTrack->GetXYZ(pos);
2262 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2263 esdTrack->GetESDpid(pid);
2264
2265 vV0FromCascade->AddDaughter(aodTrack =
2266 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2267 esdTrack->GetLabel(),
2268 p,
2269 kTRUE,
2270 pos,
2271 kFALSE,
2272 covTr,
2273 (Short_t)esdTrack->GetSign(),
2274 esdTrack->GetITSClusterMap(),
2275 pid,
2276 vV0FromCascade,
2277 kTRUE, // check if this is right
2278 kFALSE, // check if this is right
2279 AliAODTrack::kSecondary)
2280 );
2281 aodTrack->ConvertAliPIDtoAODPID();
2282 }
2283 else {
2284 cerr << "Error: event " << iEvent << " cascade " << nCascade
2285 << " track " << negFromV0 << " has already been used!" << endl;
2286 }
2287
2288 // Add the bachelor track from the cascade
2289
2290 Int_t bachelor = cascade->GetBindex();
2291
2292 if(!usedTrack[bachelor]) {
2293
2294 usedTrack[bachelor] = kTRUE;
2295
2296 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2297 esdTrack->GetPxPyPz(p);
2298 esdTrack->GetXYZ(pos);
2299 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2300 esdTrack->GetESDpid(pid);
2301
2302 vcascade->AddDaughter(aodTrack =
2303 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2304 esdTrack->GetLabel(),
2305 p,
2306 kTRUE,
2307 pos,
2308 kFALSE,
2309 covTr,
2310 (Short_t)esdTrack->GetSign(),
2311 esdTrack->GetITSClusterMap(),
2312 pid,
2313 vcascade,
2314 kTRUE, // check if this is right
2315 kFALSE, // check if this is right
2316 AliAODTrack::kSecondary)
2317 );
2318 aodTrack->ConvertAliPIDtoAODPID();
2319 }
2320 else {
2321 cerr << "Error: event " << iEvent << " cascade " << nCascade
2322 << " track " << bachelor << " has already been used!" << endl;
2323 }
2324
2325 // Add the primary track of the cascade (if any)
2326
2327 } // end of the loop on cascades
2328
2329 // V0s
2330
2331 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2332
2333 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2334
2335 AliESDv0 *v0 = esd->GetV0(nV0);
2336
2337 v0->GetXYZ(pos[0], pos[1], pos[2]);
2338 v0->GetPosCov(covVtx);
2339
2340 AliAODVertex * vV0 =
2341 new(vertices[jVertices++]) AliAODVertex(pos,
2342 covVtx,
2343 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2344 primary,
2345 AliAODVertex::kV0);
2346 primary->AddDaughter(vV0);
2347
2348 Int_t posFromV0 = v0->GetPindex();
2349 Int_t negFromV0 = v0->GetNindex();
2350
2351 // Add the positive tracks from the V0
2352
2353 if (!usedTrack[posFromV0]) {
2354
2355 usedTrack[posFromV0] = kTRUE;
2356
2357 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2358 esdTrack->GetPxPyPz(p);
2359 esdTrack->GetXYZ(pos);
2360 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2361 esdTrack->GetESDpid(pid);
2362
2363 vV0->AddDaughter(aodTrack =
2364 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2365 esdTrack->GetLabel(),
2366 p,
2367 kTRUE,
2368 pos,
2369 kFALSE,
2370 covTr,
2371 (Short_t)esdTrack->GetSign(),
2372 esdTrack->GetITSClusterMap(),
2373 pid,
2374 vV0,
2375 kTRUE, // check if this is right
2376 kFALSE, // check if this is right
2377 AliAODTrack::kSecondary)
2378 );
2379 aodTrack->ConvertAliPIDtoAODPID();
2380 }
2381 else {
2382 cerr << "Error: event " << iEvent << " V0 " << nV0
2383 << " track " << posFromV0 << " has already been used!" << endl;
2384 }
a7807689 2385
f29f1726 2386 // Add the negative tracks from the V0
2387
2388 if (!usedTrack[negFromV0]) {
2389
2390 usedTrack[negFromV0] = kTRUE;
2391
2392 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2393 esdTrack->GetPxPyPz(p);
2394 esdTrack->GetXYZ(pos);
2395 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2396 esdTrack->GetESDpid(pid);
2397
2398 vV0->AddDaughter(aodTrack =
2399 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2400 esdTrack->GetLabel(),
2401 p,
2402 kTRUE,
2403 pos,
2404 kFALSE,
2405 covTr,
2406 (Short_t)esdTrack->GetSign(),
2407 esdTrack->GetITSClusterMap(),
2408 pid,
2409 vV0,
2410 kTRUE, // check if this is right
2411 kFALSE, // check if this is right
2412 AliAODTrack::kSecondary)
2413 );
2414 aodTrack->ConvertAliPIDtoAODPID();
2415 }
2416 else {
2417 cerr << "Error: event " << iEvent << " V0 " << nV0
2418 << " track " << negFromV0 << " has already been used!" << endl;
2419 }
2420
2421 } // end of the loop on V0s
2422
2423 // Kinks: it is a big mess the access to the information in the kinks
2424 // The loop is on the tracks in order to find the mother and daugther of each kink
2425
2426
2427 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2428
2429
2430 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2431
2432 Int_t ikink = esdTrack->GetKinkIndex(0);
2433
2434 if (ikink) {
2435 // Negative kink index: mother, positive: daughter
2436
2437 // Search for the second track of the kink
2438
2439 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2440
2441 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2442
2443 Int_t jkink = esdTrack1->GetKinkIndex(0);
2444
2445 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2446
2447 // The two tracks are from the same kink
2448
2449 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2450
2451 Int_t imother = -1;
2452 Int_t idaughter = -1;
2453
2454 if (ikink<0 && jkink>0) {
2455
2456 imother = iTrack;
2457 idaughter = jTrack;
2458 }
2459 else if (ikink>0 && jkink<0) {
2460
2461 imother = jTrack;
2462 idaughter = iTrack;
2463 }
2464 else {
2465 cerr << "Error: Wrong combination of kink indexes: "
2466 << ikink << " " << jkink << endl;
2467 continue;
2468 }
2469
2470 // Add the mother track
2471
2472 AliAODTrack * mother = NULL;
2473
2474 if (!usedTrack[imother]) {
2475
2476 usedTrack[imother] = kTRUE;
2477
2478 AliESDtrack *esdTrack = esd->GetTrack(imother);
2479 esdTrack->GetPxPyPz(p);
2480 esdTrack->GetXYZ(pos);
2481 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2482 esdTrack->GetESDpid(pid);
2483
2484 mother =
2485 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2486 esdTrack->GetLabel(),
2487 p,
2488 kTRUE,
2489 pos,
2490 kFALSE,
2491 covTr,
2492 (Short_t)esdTrack->GetSign(),
2493 esdTrack->GetITSClusterMap(),
2494 pid,
2495 primary,
2496 kTRUE, // check if this is right
2497 kTRUE, // check if this is right
2498 AliAODTrack::kPrimary);
2499 primary->AddDaughter(mother);
2500 mother->ConvertAliPIDtoAODPID();
2501 }
2502 else {
2503 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2504 << " track " << imother << " has already been used!" << endl;
2505 }
2506
2507 // Add the kink vertex
2508 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2509
2510 AliAODVertex * vkink =
2511 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2512 NULL,
2513 0.,
2514 mother,
2515 AliAODVertex::kKink);
2516 // Add the daughter track
2517
2518 AliAODTrack * daughter = NULL;
2519
2520 if (!usedTrack[idaughter]) {
2521
2522 usedTrack[idaughter] = kTRUE;
2523
2524 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2525 esdTrack->GetPxPyPz(p);
2526 esdTrack->GetXYZ(pos);
2527 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2528 esdTrack->GetESDpid(pid);
2529
2530 daughter =
2531 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2532 esdTrack->GetLabel(),
2533 p,
2534 kTRUE,
2535 pos,
2536 kFALSE,
2537 covTr,
2538 (Short_t)esdTrack->GetSign(),
2539 esdTrack->GetITSClusterMap(),
2540 pid,
2541 vkink,
2542 kTRUE, // check if this is right
2543 kTRUE, // check if this is right
2544 AliAODTrack::kPrimary);
2545 vkink->AddDaughter(daughter);
2546 daughter->ConvertAliPIDtoAODPID();
2547 }
2548 else {
2549 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2550 << " track " << idaughter << " has already been used!" << endl;
2551 }
2552
2553
2554 }
2555 }
2556
2557 }
2558
2559 }
2560
2561
2562 // Tracks (primary and orphan)
2563
2564 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2565
2566
2567 if (usedTrack[nTrack]) continue;
2568
2569 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2570 esdTrack->GetPxPyPz(p);
2571 esdTrack->GetXYZ(pos);
2572 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2573 esdTrack->GetESDpid(pid);
2574
2575 Float_t impactXY, impactZ;
2576
2577 esdTrack->GetImpactParameters(impactXY,impactZ);
2578
2579 if (impactXY<3) {
2580 // track inside the beam pipe
2581
2582 primary->AddDaughter(aodTrack =
2583 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2584 esdTrack->GetLabel(),
2585 p,
2586 kTRUE,
2587 pos,
2588 kFALSE,
2589 covTr,
2590 (Short_t)esdTrack->GetSign(),
2591 esdTrack->GetITSClusterMap(),
2592 pid,
2593 primary,
2594 kTRUE, // check if this is right
2595 kTRUE, // check if this is right
2596 AliAODTrack::kPrimary)
2597 );
2598 aodTrack->ConvertAliPIDtoAODPID();
2599 }
2600 else {
2601 // outside the beam pipe: orphan track
2602 aodTrack =
2603 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2604 esdTrack->GetLabel(),
2605 p,
2606 kTRUE,
2607 pos,
2608 kFALSE,
2609 covTr,
2610 (Short_t)esdTrack->GetSign(),
2611 esdTrack->GetITSClusterMap(),
2612 pid,
2613 NULL,
2614 kFALSE, // check if this is right
2615 kFALSE, // check if this is right
2616 AliAODTrack::kOrphan);
2617 aodTrack->ConvertAliPIDtoAODPID();
2618 }
2619 } // end of loop on tracks
2620
2621 // muon tracks
2622 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2623 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2624
2625 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2626 p[0] = esdMuTrack->Px();
2627 p[1] = esdMuTrack->Py();
2628 p[2] = esdMuTrack->Pz();
2629 pos[0] = primary->GetX();
2630 pos[1] = primary->GetY();
2631 pos[2] = primary->GetZ();
2632
2633 // has to be changed once the muon pid is provided by the ESD
2634 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2635
2636 primary->AddDaughter(
2637 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2638 0, // no label provided
2639 p,
2640 kTRUE,
2641 pos,
2642 kFALSE,
2643 NULL, // no covariance matrix provided
2644 (Short_t)-99, // no charge provided
2645 0, // no ITSClusterMap
2646 pid,
2647 primary,
2648 kTRUE, // check if this is right
2649 kTRUE, // not used for vertex fit
2650 AliAODTrack::kPrimary)
2651 );
2652 }
2653
2654 // Access to the AOD container of clusters
2655 TClonesArray &clusters = *(aod->GetClusters());
2656 Int_t jClusters=0;
2657
2658 // Calo Clusters
2659 Int_t nClusters = esd->GetNumberOfCaloClusters();
2660
2661 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2662
2663 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2664
2665 Int_t id = cluster->GetID();
2666 Int_t label = -1;
53ec9628 2667 Float_t energy = cluster->E();
2668 cluster->GetPosition(posF);
f29f1726 2669 AliAODVertex *prodVertex = primary;
2670 AliAODTrack *primTrack = NULL;
2671 Char_t ttype=AliAODCluster::kUndef;
2672
2673 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2674 else if (cluster->IsEMCAL()) {
2675
2676 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2677 ttype = AliAODCluster::kEMCALPseudoCluster;
2678 else
2679 ttype = AliAODCluster::kEMCALClusterv1;
2680
2681 }
2682
2683 new(clusters[jClusters++]) AliAODCluster(id,
2684 label,
2685 energy,
2686 pos,
2687 NULL, // no covariance matrix provided
2688 NULL, // no pid for clusters provided
2689 prodVertex,
2690 primTrack,
2691 ttype);
2692
2693 } // end of loop on calo clusters
2694
2695 delete [] usedTrack;
2696 delete [] usedV0;
2697 delete [] usedKink;
2698
2699 // fill the tree for this event
2700 aodTree->Fill();
2701 } // end of event loop
2702
2703 aodTree->GetUserInfo()->Add(aod);
2704
2705 // close ESD file
2706 esdFile->Close();
2707
2708 // write the tree to the specified file
2709 aodFile = aodTree->GetCurrentFile();
2710 aodFile->cd();
2711 aodTree->Write();
2712
a7807689 2713 return;
2714}
2715
98937d93 2716void AliReconstruction::WriteAlignmentData(AliESD* esd)
2717{
2718 // Write space-points which are then used in the alignment procedures
2719 // For the moment only ITS, TRD and TPC
2720
2721 // Load TOF clusters
d528ee75 2722 if (fTracker[3]){
2723 fLoader[3]->LoadRecPoints("read");
2724 TTree* tree = fLoader[3]->TreeR();
2725 if (!tree) {
2726 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2727 return;
2728 }
2729 fTracker[3]->LoadClusters(tree);
98937d93 2730 }
98937d93 2731 Int_t ntracks = esd->GetNumberOfTracks();
2732 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2733 {
2734 AliESDtrack *track = esd->GetTrack(itrack);
2735 Int_t nsp = 0;
ef7253ac 2736 Int_t idx[200];
98937d93 2737 for (Int_t iDet = 3; iDet >= 0; iDet--)
2738 nsp += track->GetNcls(iDet);
2739 if (nsp) {
2740 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2741 track->SetTrackPointArray(sp);
2742 Int_t isptrack = 0;
2743 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2744 AliTracker *tracker = fTracker[iDet];
2745 if (!tracker) continue;
2746 Int_t nspdet = track->GetNcls(iDet);
98937d93 2747 if (nspdet <= 0) continue;
2748 track->GetClusters(iDet,idx);
2749 AliTrackPoint p;
2750 Int_t isp = 0;
2751 Int_t isp2 = 0;
2752 while (isp < nspdet) {
2753 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 2754 const Int_t kNTPCmax = 159;
2755 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2756 if (!isvalid) continue;
2757 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2758 }
98937d93 2759 }
2760 }
2761 }
d528ee75 2762 if (fTracker[3]){
2763 fTracker[3]->UnloadClusters();
2764 fLoader[3]->UnloadRecPoints();
2765 }
98937d93 2766}
2e3550da 2767
2768//_____________________________________________________________________________
2769void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2770{
2771 // The method reads the raw-data error log
2772 // accumulated within the rawReader.
2773 // It extracts the raw-data errors related to
2774 // the current event and stores them into
2775 // a TClonesArray inside the esd object.
2776
2777 if (!fRawReader) return;
2778
2779 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2780
2781 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2782 if (!log) continue;
2783 if (iEvent != log->GetEventNumber()) continue;
2784
2785 esd->AddRawDataErrorLog(log);
2786 }
2787
2788}
46698ae4 2789
2790TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
b545009a 2791 // Dump a file content into a char in TNamed
46698ae4 2792 ifstream in;
2793 in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2794 Int_t kBytes = (Int_t)in.tellg();
2795 printf("Size: %d \n",kBytes);
2796 TNamed *fn = 0;
2797 if(in.good()){
2798 char* memblock = new char [kBytes];
2799 in.seekg (0, ios::beg);
2800 in.read (memblock, kBytes);
2801 in.close();
2802 TString fData(memblock,kBytes);
2803 fn = new TNamed(fName,fData);
2804 printf("fData Size: %d \n",fData.Sizeof());
2805 printf("fName Size: %d \n",fName.Sizeof());
2806 printf("fn Size: %d \n",fn->Sizeof());
2807 delete[] memblock;
2808 }
2809 else{
2810 AliInfo(Form("Could not Open %s\n",fPath.Data()));
2811 }
2812
2813 return fn;
2814}
2815
2816void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
46698ae4 2817 // This is not really needed in AliReconstruction at the moment
2818 // but can serve as a template
2819
2820 TList *fList = fTree->GetUserInfo();
2821 TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2822 printf("fn Size: %d \n",fn->Sizeof());
2823
2824 TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2825 const char* cdata = fn->GetTitle();
2826 printf("fTmp Size %d\n",fTmp.Sizeof());
2827
2828 int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2829 printf("calculated size %d\n",size);
2830 ofstream out(fName.Data(),ios::out | ios::binary);
2831 out.write(cdata,size);
2832 out.close();
2833
2834}
2835