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