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