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