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