]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliReconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
corrected bug in GetValue(UInt_t timeSec) (Haavard)
[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"
596a855f 135#include "AliESD.h"
1d99986f 136#include "AliESDfriend.h"
2257f27e 137#include "AliESDVertex.h"
32e449be 138#include "AliMultiplicity.h"
c84a5e9e 139#include "AliTracker.h"
2257f27e 140#include "AliVertexer.h"
c5e3e5d1 141#include "AliVertexerTracks.h"
5e4ff34d 142#include "AliV0vertexer.h"
143#include "AliCascadeVertexer.h"
596a855f 144#include "AliHeader.h"
145#include "AliGenEventHeader.h"
b26c3770 146#include "AliPID.h"
596a855f 147#include "AliESDpid.h"
ff8bb5ae 148#include "AliESDtrack.h"
f3a97c86 149
150#include "AliRunTag.h"
f3a97c86 151#include "AliDetectorTag.h"
152#include "AliEventTag.h"
153
25be1e5c 154#include "AliGeomManager.h"
98937d93 155#include "AliTrackPointArray.h"
b0314964 156#include "AliCDBManager.h"
6bae477a 157#include "AliCDBEntry.h"
158#include "AliAlignObj.h"
f3a97c86 159
b647652d 160#include "AliCentralTrigger.h"
161#include "AliCTPRawStream.h"
162
f29f1726 163#include "AliAODEvent.h"
164#include "AliAODHeader.h"
165#include "AliAODTrack.h"
166#include "AliAODVertex.h"
167#include "AliAODCluster.h"
168
169
596a855f 170ClassImp(AliReconstruction)
171
172
c757bafd 173//_____________________________________________________________________________
b384f8a4 174const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
c757bafd 175
596a855f 176//_____________________________________________________________________________
024cf675 177AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 178 const char* name, const char* title) :
179 TNamed(name, title),
180
c84a5e9e 181 fUniformField(kTRUE),
2257f27e 182 fRunVertexFinder(kTRUE),
1f46a9ae 183 fRunHLTTracking(kFALSE),
e66fbafb 184 fRunMuonTracking(kFALSE),
1d99986f 185 fStopOnError(kFALSE),
186 fWriteAlignmentData(kFALSE),
187 fWriteESDfriend(kFALSE),
a7807689 188 fWriteAOD(kFALSE),
b647652d 189 fFillTriggerESD(kTRUE),
1d99986f 190
191 fRunLocalReconstruction("ALL"),
b8cd5251 192 fRunTracking("ALL"),
e583c30d 193 fFillESD("ALL"),
194 fGAliceFileName(gAliceFilename),
b649205a 195 fInput(""),
35042093 196 fEquipIdMap(""),
b26c3770 197 fFirstEvent(0),
198 fLastEvent(-1),
973388c2 199 fNumberOfEventsPerFile(1),
24f7a148 200 fCheckPointLevel(0),
b8cd5251 201 fOptions(),
6bae477a 202 fLoadAlignFromCDB(kTRUE),
203 fLoadAlignData("ALL"),
46698ae4 204 fESDPar(""),
e583c30d 205
206 fRunLoader(NULL),
b649205a 207 fRawReader(NULL),
b8cd5251 208
98937d93 209 fVertexer(NULL),
9178838a 210 fDiamondProfile(NULL),
98937d93 211
6bae477a 212 fAlignObjArray(NULL),
ec92bee0 213 fCDBUri(cdbUri),
214 fSpecCDBUri()
596a855f 215{
216// create reconstruction object with default parameters
b8cd5251 217
218 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
219 fReconstructor[iDet] = NULL;
220 fLoader[iDet] = NULL;
221 fTracker[iDet] = NULL;
222 }
e47c4c2e 223 AliPID pid;
596a855f 224}
225
226//_____________________________________________________________________________
227AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 228 TNamed(rec),
229
c84a5e9e 230 fUniformField(rec.fUniformField),
2257f27e 231 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 232 fRunHLTTracking(rec.fRunHLTTracking),
e66fbafb 233 fRunMuonTracking(rec.fRunMuonTracking),
1d99986f 234 fStopOnError(rec.fStopOnError),
235 fWriteAlignmentData(rec.fWriteAlignmentData),
236 fWriteESDfriend(rec.fWriteESDfriend),
a7807689 237 fWriteAOD(rec.fWriteAOD),
b647652d 238 fFillTriggerESD(rec.fFillTriggerESD),
1d99986f 239
240 fRunLocalReconstruction(rec.fRunLocalReconstruction),
e583c30d 241 fRunTracking(rec.fRunTracking),
242 fFillESD(rec.fFillESD),
243 fGAliceFileName(rec.fGAliceFileName),
b649205a 244 fInput(rec.fInput),
35042093 245 fEquipIdMap(rec.fEquipIdMap),
b26c3770 246 fFirstEvent(rec.fFirstEvent),
247 fLastEvent(rec.fLastEvent),
973388c2 248 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
24f7a148 249 fCheckPointLevel(0),
b8cd5251 250 fOptions(),
6bae477a 251 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
252 fLoadAlignData(rec.fLoadAlignData),
46698ae4 253 fESDPar(rec.fESDPar),
e583c30d 254
255 fRunLoader(NULL),
b649205a 256 fRawReader(NULL),
b8cd5251 257
98937d93 258 fVertexer(NULL),
9178838a 259 fDiamondProfile(NULL),
98937d93 260
6bae477a 261 fAlignObjArray(rec.fAlignObjArray),
ec92bee0 262 fCDBUri(rec.fCDBUri),
263 fSpecCDBUri()
596a855f 264{
265// copy constructor
266
ec92bee0 267 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
efd2085e 268 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
269 }
b8cd5251 270 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
271 fReconstructor[iDet] = NULL;
272 fLoader[iDet] = NULL;
273 fTracker[iDet] = NULL;
274 }
ec92bee0 275 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
276 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
277 }
596a855f 278}
279
280//_____________________________________________________________________________
281AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
282{
283// assignment operator
284
285 this->~AliReconstruction();
286 new(this) AliReconstruction(rec);
287 return *this;
288}
289
290//_____________________________________________________________________________
291AliReconstruction::~AliReconstruction()
292{
293// clean up
294
e583c30d 295 CleanUp();
efd2085e 296 fOptions.Delete();
ec92bee0 297 fSpecCDBUri.Delete();
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
6bae477a 428//_____________________________________________________________________________
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
46698ae4 562 AliESD* esd = new AliESD(); AliESD* hltesd = new AliESD();
b26c3770 563 TFile* fileOld = NULL;
1f46a9ae 564 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 565 if (!gSystem->AccessPathName("AliESDs.root")){
566 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
567 fileOld = TFile::Open("AliESDs.old.root");
568 if (fileOld && fileOld->IsOpen()) {
569 treeOld = (TTree*) fileOld->Get("esdTree");
46698ae4 570 if (treeOld)esd->ReadFromTree(treeOld);
1f46a9ae 571 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
46698ae4 572 if (hlttreeOld) hltesd->ReadFromTree(hlttreeOld);
b26c3770 573 }
574 }
575
36711aa4 576 // create the ESD output file and tree
596a855f 577 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
46698ae4 578 file->SetCompressionLevel(2);
596a855f 579 if (!file->IsOpen()) {
815c2b38 580 AliError("opening AliESDs.root failed");
b26c3770 581 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 582 }
46698ae4 583
36711aa4 584 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
46698ae4 585 esd = new AliESD();
586 esd->CreateStdContent();
587 esd->WriteToTree(tree);
588
1f46a9ae 589 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
46698ae4 590 hltesd = new AliESD();
591 hltesd->CreateStdContent();
592 hltesd->WriteToTree(hlttree);
593
594 /* CKB Why?
1f46a9ae 595 delete esd; delete hltesd;
596 esd = NULL; hltesd = NULL;
46698ae4 597 */
500d54ab 598 // create the branch with ESD additions
46698ae4 599 AliESDfriend *esdf = 0;
1d99986f 600 if (fWriteESDfriend) {
46698ae4 601 esdf = new AliESDfriend();
602 TBranch *br=tree->Branch("ESDfriend.","AliESDfriend", &esdf);
603 br->SetFile("AliESDfriends.root");
604 esd->AddObject(esdf);
1d99986f 605 }
46698ae4 606
17c86e90 607 // Get the diamond profile from OCDB
608 AliCDBEntry* entry = AliCDBManager::Instance()
609 ->Get("GRP/Calib/MeanVertex");
610
611 if(entry) {
612 fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
613 } else {
614 AliError("No diamond profile found in OCDB!");
615 }
616
20e5681c 617 AliVertexerTracks tVertexer(AliTracker::GetBz());
9178838a 618 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
c5e3e5d1 619
596a855f 620 // loop over events
b649205a 621 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 622
596a855f 623 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 624 if (fRawReader) fRawReader->NextEvent();
4a33489c 625 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
b26c3770 626 // copy old ESD to the new one
627 if (treeOld) {
46698ae4 628 esd->ReadFromTree(treeOld);
b26c3770 629 treeOld->GetEntry(iEvent);
630 }
631 tree->Fill();
1f46a9ae 632 if (hlttreeOld) {
46698ae4 633 esd->ReadFromTree(hlttreeOld);
1f46a9ae 634 hlttreeOld->GetEntry(iEvent);
635 }
636 hlttree->Fill();
b26c3770 637 continue;
638 }
46698ae4 639
815c2b38 640 AliInfo(Form("processing event %d", iEvent));
596a855f 641 fRunLoader->GetEvent(iEvent);
24f7a148 642
bb0901a4 643 char aFileName[256];
644 sprintf(aFileName, "ESD_%d.%d_final.root",
f08fc9f5 645 fRunLoader->GetHeader()->GetRun(),
646 fRunLoader->GetHeader()->GetEventNrInRun());
bb0901a4 647 if (!gSystem->AccessPathName(aFileName)) continue;
24f7a148 648
b26c3770 649 // local reconstruction
650 if (!fRunLocalReconstruction.IsNull()) {
651 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
652 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
653 }
654 }
655
46698ae4 656
f08fc9f5 657 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 658 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
31fd97b2 659 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
660 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
46698ae4 661
d6ee376f 662 // Set magnetic field from the tracker
663 esd->SetMagneticField(AliTracker::GetBz());
664 hltesd->SetMagneticField(AliTracker::GetBz());
596a855f 665
46698ae4 666
667
2e3550da 668 // Fill raw-data error log into the ESD
669 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
670
2257f27e 671 // vertex finder
672 if (fRunVertexFinder) {
673 if (!ReadESD(esd, "vertex")) {
674 if (!RunVertexFinder(esd)) {
b26c3770 675 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 676 }
677 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
678 }
679 }
680
1f46a9ae 681 // HLT tracking
682 if (!fRunTracking.IsNull()) {
683 if (fRunHLTTracking) {
684 hltesd->SetVertex(esd->GetVertex());
685 if (!RunHLTTracking(hltesd)) {
686 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
687 }
688 }
689 }
690
e66fbafb 691 // Muon tracking
b8cd5251 692 if (!fRunTracking.IsNull()) {
e66fbafb 693 if (fRunMuonTracking) {
761350a6 694 if (!RunMuonTracking(esd)) {
b26c3770 695 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 696 }
596a855f 697 }
698 }
699
e66fbafb 700 // barrel tracking
701 if (!fRunTracking.IsNull()) {
21c573b7 702 if (!ReadESD(esd, "tracking")) {
703 if (!RunTracking(esd)) {
704 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
e66fbafb 705 }
21c573b7 706 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
e66fbafb 707 }
708 }
21c573b7 709
596a855f 710 // fill ESD
711 if (!fFillESD.IsNull()) {
712 if (!FillESD(esd, fFillESD)) {
b26c3770 713 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 714 }
715 }
001397cd 716 // fill Event header information from the RawEventHeader
717 if (fRawReader){FillRawEventHeaderESD(esd);}
596a855f 718
719 // combined PID
720 AliESDpid::MakePID(esd);
24f7a148 721 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 722
b647652d 723 if (fFillTriggerESD) {
724 if (!ReadESD(esd, "trigger")) {
725 if (!FillTriggerESD(esd)) {
726 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
727 }
728 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
729 }
730 }
731
a03cd2e0 732 //Try to improve the reconstructed primary vertex position using the tracks
c060d7fe 733 AliESDVertex *pvtx=0;
734 Bool_t dovertex=kTRUE;
735 TObject* obj = fOptions.FindObject("ITS");
736 if (obj) {
737 TString optITS = obj->GetTitle();
738 if (optITS.Contains("cosmics") || optITS.Contains("COSMICS"))
739 dovertex=kFALSE;
740 }
741 if(dovertex) pvtx=tVertexer.FindPrimaryVertex(esd);
17c86e90 742 if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
743
a03cd2e0 744 if (pvtx)
745 if (pvtx->GetStatus()) {
746 // Store the improved primary vertex
747 esd->SetPrimaryVertex(pvtx);
748 // Propagate the tracks to the DCA to the improved primary vertex
749 Double_t somethingbig = 777.;
750 Double_t bz = esd->GetMagneticField();
751 Int_t nt=esd->GetNumberOfTracks();
752 while (nt--) {
753 AliESDtrack *t = esd->GetTrack(nt);
754 t->RelateToVertex(pvtx, bz, somethingbig);
755 }
756 }
c5e3e5d1 757
5e4ff34d 758 {
759 // V0 finding
760 AliV0vertexer vtxer;
761 vtxer.Tracks2V0vertices(esd);
762
763 // Cascade finding
764 AliCascadeVertexer cvtxer;
765 cvtxer.V0sTracks2CascadeVertices(esd);
766 }
767
596a855f 768 // write ESD
1d99986f 769 if (fWriteESDfriend) {
46698ae4 770 new (esdf) AliESDfriend(); // Reset...
771 esd->GetESDfriend(esdf);
1d99986f 772 }
500d54ab 773 tree->Fill();
774
775 // write HLT ESD
776 hlttree->Fill();
1d99986f 777
f3a97c86 778 if (fCheckPointLevel > 0) WriteESD(esd, "final");
46698ae4 779 esd->Reset();
780 hltesd->Reset();
781 // esdf->Reset();
782 // delete esdf; esdf = 0;
783 }
784
785
786 tree->GetUserInfo()->Add(esd);
787 hlttree->GetUserInfo()->Add(hltesd);
788
789
790
791 if(fESDPar.Contains("ESD.par")){
792 AliInfo("Attaching ESD.par to Tree");
793 TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
794 tree->GetUserInfo()->Add(fn);
596a855f 795 }
796
46698ae4 797
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
b26c3770 875//_____________________________________________________________________________
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
596a855f 938//_____________________________________________________________________________
2257f27e 939Bool_t AliReconstruction::RunVertexFinder(AliESD*& 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
1f46a9ae 994//_____________________________________________________________________________
995Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
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
e66fbafb 1050//_____________________________________________________________________________
761350a6 1051Bool_t AliReconstruction::RunMuonTracking(AliESD*& 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
2257f27e 1110//_____________________________________________________________________________
1111Bool_t AliReconstruction::RunTracking(AliESD*& esd)
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//_____________________________________________________________________________
24f7a148 1235Bool_t AliReconstruction::FillESD(AliESD*& 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//_____________________________________________________________________________
1298Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
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//_____________________________________________________________________________
1344Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
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
f08fc9f5 1405//_____________________________________________________________________________
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
c757bafd 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
2257f27e 1573//_____________________________________________________________________________
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
24f7a148 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
e583c30d 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//_____________________________________________________________________________
1658Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1659{
1660// read the ESD event from a file
1661
1662 if (!esd) return kFALSE;
1663 char fileName[256];
1664 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1665 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1666 if (gSystem->AccessPathName(fileName)) return kFALSE;
1667
f3a97c86 1668 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1669 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1670 TFile* file = TFile::Open(fileName);
1671 if (!file || !file->IsOpen()) {
815c2b38 1672 AliError(Form("opening %s failed", fileName));
24f7a148 1673 delete file;
1674 return kFALSE;
1675 }
1676
1677 gROOT->cd();
1678 delete esd;
1679 esd = (AliESD*) file->Get("ESD");
1680 file->Close();
1681 delete file;
1682 return kTRUE;
1683}
1684
1685//_____________________________________________________________________________
1686void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1687{
1688// write the ESD event to a file
1689
1690 if (!esd) return;
1691 char fileName[256];
1692 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1693 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1694
815c2b38 1695 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1696 TFile* file = TFile::Open(fileName, "recreate");
1697 if (!file || !file->IsOpen()) {
815c2b38 1698 AliError(Form("opening %s failed", fileName));
24f7a148 1699 } else {
1700 esd->Write("ESD");
1701 file->Close();
1702 }
1703 delete file;
1704}
f3a97c86 1705
1706
1707
1708
1709//_____________________________________________________________________________
1710void AliReconstruction::CreateTag(TFile* file)
1711{
1ac7d449 1712 //GRP
1713 Float_t lhcLuminosity = 0.0;
1714 TString lhcState = "test";
1715 UInt_t detectorMask = 0;
1716
2bdb9d38 1717 /////////////
1718 //muon code//
1719 ////////////
56982dd1 1720 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1721 //Variables
56982dd1 1722 Double_t fX,fY,fZ ;
1723 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1724 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1725 Int_t fCharge;
1726 TLorentzVector fEPvector;
1727
1728 Float_t fZVertexCut = 10.0;
1729 Float_t fRhoVertexCut = 2.0;
1730
1731 Float_t fLowPtCut = 1.0;
1732 Float_t fHighPtCut = 3.0;
1733 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1734 ////////////
1735
1736 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1737
089bf903 1738 // Creates the tags for all the events in a given ESD file
4302e20f 1739 Int_t ntrack;
089bf903 1740 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1741 Int_t nPos, nNeg, nNeutr;
1742 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1743 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1744 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1745 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1746 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1747 Int_t fVertexflag;
1387d165 1748 Int_t iRunNumber = 0;
547d75a6 1749 TString fVertexName("default");
4302e20f 1750
f3a97c86 1751 AliRunTag *tag = new AliRunTag();
f3a97c86 1752 AliEventTag *evTag = new AliEventTag();
1753 TTree ttag("T","A Tree with event tags");
2bdb9d38 1754 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1755 btag->SetCompressionLevel(9);
2bdb9d38 1756
f3a97c86 1757 AliInfo(Form("Creating the tags......."));
1758
1759 if (!file || !file->IsOpen()) {
1760 AliError(Form("opening failed"));
1761 delete file;
1762 return ;
2bdb9d38 1763 }
d71aa170 1764 Int_t lastEvent = 0;
46698ae4 1765 TTree *b = (TTree*) file->Get("esdTree");
1766 AliESD *esd = new AliESD();
1767 esd->ReadFromTree(b);
1387d165 1768
d71aa170 1769 b->GetEntry(fFirstEvent);
1387d165 1770 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1771
bb0901a4 1772 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1773 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1774 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 1775 ntrack = 0;
1776 nPos = 0;
1777 nNeg = 0;
1778 nNeutr =0;
1779 nK0s = 0;
1780 nNeutrons = 0;
1781 nPi0s = 0;
1782 nGamas = 0;
1783 nProtons = 0;
1784 nKaons = 0;
1785 nPions = 0;
1786 nMuons = 0;
1787 nElectrons = 0;
1788 nCh1GeV = 0;
1789 nCh3GeV = 0;
1790 nCh10GeV = 0;
1791 nMu1GeV = 0;
1792 nMu3GeV = 0;
1793 nMu10GeV = 0;
1794 nEl1GeV = 0;
1795 nEl3GeV = 0;
1796 nEl10GeV = 0;
1797 maxPt = .0;
1798 meanPt = .0;
1799 totalP = .0;
547d75a6 1800 fVertexflag = 0;
d1a50cb5 1801
2bdb9d38 1802 b->GetEntry(iEventNumber);
1387d165 1803 iRunNumber = esd->GetRunNumber();
1804 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1805 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1806 if (!vertexIn) AliError("ESD has not defined vertex.");
1807 if (vertexIn) fVertexName = vertexIn->GetName();
1808 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1809 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1810 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1811 UInt_t status = esdTrack->GetStatus();
f3a97c86 1812
2bdb9d38 1813 //select only tracks with ITS refit
1814 if ((status&AliESDtrack::kITSrefit)==0) continue;
1815 //select only tracks with TPC refit
1816 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1817
2bdb9d38 1818 //select only tracks with the "combined PID"
1819 if ((status&AliESDtrack::kESDpid)==0) continue;
1820 Double_t p[3];
1821 esdTrack->GetPxPyPz(p);
1822 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1823 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1824 totalP += momentum;
1825 meanPt += fPt;
1826 if(fPt > maxPt) maxPt = fPt;
4302e20f 1827
2bdb9d38 1828 if(esdTrack->GetSign() > 0) {
1829 nPos++;
56982dd1 1830 if(fPt > fLowPtCut) nCh1GeV++;
1831 if(fPt > fHighPtCut) nCh3GeV++;
1832 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1833 }
1834 if(esdTrack->GetSign() < 0) {
1835 nNeg++;
56982dd1 1836 if(fPt > fLowPtCut) nCh1GeV++;
1837 if(fPt > fHighPtCut) nCh3GeV++;
1838 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1839 }
1840 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1841
2bdb9d38 1842 //PID
1843 Double_t prob[5];
1844 esdTrack->GetESDpid(prob);
4302e20f 1845
2bdb9d38 1846 Double_t rcc = 0.0;
1847 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1848 if(rcc == 0.0) continue;
1849 //Bayes' formula
1850 Double_t w[5];
1851 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1852
2bdb9d38 1853 //protons
1854 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1855 //kaons
1856 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1857 //pions
1858 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1859 //electrons
1860 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1861 nElectrons++;
56982dd1 1862 if(fPt > fLowPtCut) nEl1GeV++;
1863 if(fPt > fHighPtCut) nEl3GeV++;
1864 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1865 }
1866 ntrack++;
1867 }//track loop
1868
1869 /////////////
1870 //muon code//
1871 ////////////
1872 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1873 // loop over all reconstructed tracks (also first track of combination)
1874 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1875 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1876 if (muonTrack == 0x0) continue;
4302e20f 1877
2bdb9d38 1878 // Coordinates at vertex
56982dd1 1879 fZ = muonTrack->GetZ();
1880 fY = muonTrack->GetBendingCoor();
1881 fX = muonTrack->GetNonBendingCoor();
4302e20f 1882
56982dd1 1883 fThetaX = muonTrack->GetThetaX();
1884 fThetaY = muonTrack->GetThetaY();
4302e20f 1885
56982dd1 1886 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1887 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1888 fPxRec = fPzRec * TMath::Tan(fThetaX);
1889 fPyRec = fPzRec * TMath::Tan(fThetaY);
1890 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1891
2bdb9d38 1892 //ChiSquare of the track if needed
56982dd1 1893 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1894 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1895 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1896
2bdb9d38 1897 // total number of muons inside a vertex cut
56982dd1 1898 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1899 nMuons++;
56982dd1 1900 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1901 nMu1GeV++;
56982dd1 1902 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1903 nMu3GeV++;
56982dd1 1904 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1905 nMu10GeV++;
1906 }
1907 }
1908 }
1909 }
1910 }//muon track loop
1911
1912 // Fill the event tags
1913 if(ntrack != 0)
1914 meanPt = meanPt/ntrack;
1915
1916 evTag->SetEventId(iEventNumber+1);
547d75a6 1917 if (vertexIn) {
1918 evTag->SetVertexX(vertexIn->GetXv());
1919 evTag->SetVertexY(vertexIn->GetYv());
1920 evTag->SetVertexZ(vertexIn->GetZv());
1921 evTag->SetVertexZError(vertexIn->GetZRes());
1922 }
d1a50cb5 1923 evTag->SetVertexFlag(fVertexflag);
1924
2bdb9d38 1925 evTag->SetT0VertexZ(esd->GetT0zVertex());
1926
8bd8ac26 1927 evTag->SetTriggerMask(esd->GetTriggerMask());
1928 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1929
32a5cab4 1930 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1931 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1932 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1933 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1934 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1935 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1936
1937
1938 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1939 evTag->SetNumOfPosTracks(nPos);
1940 evTag->SetNumOfNegTracks(nNeg);
1941 evTag->SetNumOfNeutrTracks(nNeutr);
1942
1943 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1944 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1945 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1946 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1947
1948 evTag->SetNumOfProtons(nProtons);
1949 evTag->SetNumOfKaons(nKaons);
1950 evTag->SetNumOfPions(nPions);
1951 evTag->SetNumOfMuons(nMuons);
1952 evTag->SetNumOfElectrons(nElectrons);
1953 evTag->SetNumOfPhotons(nGamas);
1954 evTag->SetNumOfPi0s(nPi0s);
1955 evTag->SetNumOfNeutrons(nNeutrons);
1956 evTag->SetNumOfKaon0s(nK0s);
1957
1958 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1959 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1960 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1961 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1962 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1963 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1964 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1965 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1966 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1967
85c60a8e 1968 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1969 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1970
1971 evTag->SetTotalMomentum(totalP);
1972 evTag->SetMeanPt(meanPt);
1973 evTag->SetMaxPt(maxPt);
1974
1ac7d449 1975 tag->SetLHCTag(lhcLuminosity,lhcState);
1976 tag->SetDetectorTag(detectorMask);
1977
1387d165 1978 tag->SetRunId(iInitRunNumber);
2bdb9d38 1979 tag->AddEventTag(*evTag);
1980 }
bb0901a4 1981 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 1982 else lastEvent = fLastEvent;
f3a97c86 1983
1984 ttag.Fill();
1985 tag->Clear();
1986
1987 char fileName[256];
1988 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 1989 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 1990 AliInfo(Form("writing tags to file %s", fileName));
1991 AliDebug(1, Form("writing tags to file %s", fileName));
1992
1993 TFile* ftag = TFile::Open(fileName, "recreate");
1994 ftag->cd();
1995 ttag.Write();
1996 ftag->Close();
1997 file->cd();
1998 delete tag;
f3a97c86 1999 delete evTag;
2000}
2001
a7807689 2002//_____________________________________________________________________________
f29f1726 2003void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
a7807689 2004{
f29f1726 2005 // write all files from the given esd file to an aod file
2006
2007 // create an AliAOD object
2008 AliAODEvent *aod = new AliAODEvent();
2009 aod->CreateStdContent();
2010
2011 // go to the file
2012 aodFile->cd();
2013
2014 // create the tree
2015 TTree *aodTree = new TTree("AOD", "AliAOD tree");
2016 aodTree->Branch(aod->GetList());
2017
2018 // connect to ESD
2019 TTree *t = (TTree*) esdFile->Get("esdTree");
53ec9628 2020 AliESD *esd = new AliESD();
2021 esd->ReadFromTree(t);
f29f1726 2022
53ec9628 2023 Int_t nEvents = t->GetEntries();
f29f1726 2024
2025 // set arrays and pointers
2026 Float_t posF[3];
2027 Double_t pos[3];
2028 Double_t p[3];
2029 Double_t covVtx[6];
2030 Double_t covTr[21];
2031 Double_t pid[10];
2032
2033 // loop over events and fill them
2034 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
53ec9628 2035 t->GetEntry(iEvent);
f29f1726 2036
2037 // Multiplicity information needed by the header (to be revised!)
2038 Int_t nTracks = esd->GetNumberOfTracks();
2039 Int_t nPosTracks = 0;
2040 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2041 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2042
a1d4139d 2043 // Update the header
2044 AliAODHeader* header = aod->GetHeader();
2045
ade23daf 2046 header->SetRunNumber (esd->GetRunNumber() );
2047 header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
2048 header->SetOrbitNumber (esd->GetOrbitNumber() );
2049 header->SetPeriodNumber (esd->GetPeriodNumber() );
2050 header->SetTriggerMask (esd->GetTriggerMask() );
2051 header->SetTriggerCluster (esd->GetTriggerCluster() );
2052 header->SetEventType (esd->GetEventType() );
2053 header->SetMagneticField (esd->GetMagneticField() );
2054 header->SetZDCN1Energy (esd->GetZDCN1Energy() );
2055 header->SetZDCP1Energy (esd->GetZDCP1Energy() );
2056 header->SetZDCN2Energy (esd->GetZDCN2Energy() );
2057 header->SetZDCP2Energy (esd->GetZDCP2Energy() );
2058 header->SetZDCEMEnergy (esd->GetZDCEMEnergy() );
a1d4139d 2059 header->SetRefMultiplicity (nTracks);
2060 header->SetRefMultiplicityPos(nPosTracks);
2061 header->SetRefMultiplicityNeg(nTracks - nPosTracks);
2062 header->SetMuonMagFieldScale(-999.); // FIXME
2063 header->SetCentrality(-999.); // FIXME
2064//
2065//
f29f1726 2066
2067 Int_t nV0s = esd->GetNumberOfV0s();
2068 Int_t nCascades = esd->GetNumberOfCascades();
2069 Int_t nKinks = esd->GetNumberOfKinks();
2070 Int_t nVertices = nV0s + nCascades + nKinks;
2071
2072 aod->ResetStd(nTracks, nVertices);
2073 AliAODTrack *aodTrack;
2074
2075
2076 // Array to take into account the tracks already added to the AOD
2077 Bool_t * usedTrack = NULL;
2078 if (nTracks>0) {
2079 usedTrack = new Bool_t[nTracks];
2080 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2081 }
2082 // Array to take into account the V0s already added to the AOD
2083 Bool_t * usedV0 = NULL;
2084 if (nV0s>0) {
2085 usedV0 = new Bool_t[nV0s];
2086 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2087 }
2088 // Array to take into account the kinks already added to the AOD
2089 Bool_t * usedKink = NULL;
2090 if (nKinks>0) {
2091 usedKink = new Bool_t[nKinks];
2092 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2093 }
2094
2095 // Access to the AOD container of vertices
2096 TClonesArray &vertices = *(aod->GetVertices());
2097 Int_t jVertices=0;
2098
2099 // Access to the AOD container of tracks
2100 TClonesArray &tracks = *(aod->GetTracks());
2101 Int_t jTracks=0;
2102
2103 // Add primary vertex. The primary tracks will be defined
2104 // after the loops on the composite objects (V0, cascades, kinks)
2105 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2106
2107 vtx->GetXYZ(pos); // position
2108 vtx->GetCovMatrix(covVtx); //covariance matrix
2109
2110 AliAODVertex * primary = new(vertices[jVertices++])
2111 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2112
2113 // Create vertices starting from the most complex objects
2114
2115 // Cascades
2116 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2117 AliESDcascade *cascade = esd->GetCascade(nCascade);
2118
2119 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2120 cascade->GetPosCovXi(covVtx);
2121
2122 // Add the cascade vertex
2123 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2124 covVtx,
2125 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2126 primary,
2127 AliAODVertex::kCascade);
2128
2129 primary->AddDaughter(vcascade);
2130
2131 // Add the V0 from the cascade. The ESD class have to be optimized...
2132 // Now we have to search for the corresponding Vo in the list of V0s
2133 // using the indeces of the positive and negative tracks
2134
2135 Int_t posFromV0 = cascade->GetPindex();
2136 Int_t negFromV0 = cascade->GetNindex();
2137
2138
2139 AliESDv0 * v0 = 0x0;
2140 Int_t indV0 = -1;
2141
2142 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2143
2144 v0 = esd->GetV0(iV0);
2145 Int_t posV0 = v0->GetPindex();
2146 Int_t negV0 = v0->GetNindex();
2147
2148 if (posV0==posFromV0 && negV0==negFromV0) {
2149 indV0 = iV0;
2150 break;
2151 }
2152 }
2153
2154 AliAODVertex * vV0FromCascade = 0x0;
2155
2156 if (indV0>-1 && !usedV0[indV0] ) {
2157
2158 // the V0 exists in the array of V0s and is not used
2159
2160 usedV0[indV0] = kTRUE;
2161
2162 v0->GetXYZ(pos[0], pos[1], pos[2]);
2163 v0->GetPosCov(covVtx);
2164
2165 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2166 covVtx,
2167 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2168 vcascade,
2169 AliAODVertex::kV0);
2170 } else {
2171
2172 // the V0 doesn't exist in the array of V0s or was used
2173 cerr << "Error: event " << iEvent << " cascade " << nCascade
2174 << " The V0 " << indV0
2175 << " doesn't exist in the array of V0s or was used!" << endl;
2176
2177 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2178 cascade->GetPosCov(covVtx);
2179
2180 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2181 covVtx,
2182 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2183 vcascade,
2184 AliAODVertex::kV0);
2185 vcascade->AddDaughter(vV0FromCascade);
2186 }
2187
2188 // Add the positive tracks from the V0
2189
2190 if (! usedTrack[posFromV0]) {
2191
2192 usedTrack[posFromV0] = kTRUE;
2193
2194 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2195 esdTrack->GetPxPyPz(p);
2196 esdTrack->GetXYZ(pos);
2197 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2198 esdTrack->GetESDpid(pid);
2199
2200 vV0FromCascade->AddDaughter(aodTrack =
2201 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2202 esdTrack->GetLabel(),
2203 p,
2204 kTRUE,
2205 pos,
2206 kFALSE,
2207 covTr,
2208 (Short_t)esdTrack->GetSign(),
2209 esdTrack->GetITSClusterMap(),
2210 pid,
2211 vV0FromCascade,
2212 kTRUE, // check if this is right
2213 kFALSE, // check if this is right
2214 AliAODTrack::kSecondary)
2215 );
2216 aodTrack->ConvertAliPIDtoAODPID();
2217 }
2218 else {
2219 cerr << "Error: event " << iEvent << " cascade " << nCascade
2220 << " track " << posFromV0 << " has already been used!" << endl;
2221 }
2222
2223 // Add the negative tracks from the V0
2224
2225 if (!usedTrack[negFromV0]) {
2226
2227 usedTrack[negFromV0] = kTRUE;
2228
2229 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2230 esdTrack->GetPxPyPz(p);
2231 esdTrack->GetXYZ(pos);
2232 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2233 esdTrack->GetESDpid(pid);
2234
2235 vV0FromCascade->AddDaughter(aodTrack =
2236 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2237 esdTrack->GetLabel(),
2238 p,
2239 kTRUE,
2240 pos,
2241 kFALSE,
2242 covTr,
2243 (Short_t)esdTrack->GetSign(),
2244 esdTrack->GetITSClusterMap(),
2245 pid,
2246 vV0FromCascade,
2247 kTRUE, // check if this is right
2248 kFALSE, // check if this is right
2249 AliAODTrack::kSecondary)
2250 );
2251 aodTrack->ConvertAliPIDtoAODPID();
2252 }
2253 else {
2254 cerr << "Error: event " << iEvent << " cascade " << nCascade
2255 << " track " << negFromV0 << " has already been used!" << endl;
2256 }
2257
2258 // Add the bachelor track from the cascade
2259
2260 Int_t bachelor = cascade->GetBindex();
2261
2262 if(!usedTrack[bachelor]) {
2263
2264 usedTrack[bachelor] = kTRUE;
2265
2266 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2267 esdTrack->GetPxPyPz(p);
2268 esdTrack->GetXYZ(pos);
2269 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2270 esdTrack->GetESDpid(pid);
2271
2272 vcascade->AddDaughter(aodTrack =
2273 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2274 esdTrack->GetLabel(),
2275 p,
2276 kTRUE,
2277 pos,
2278 kFALSE,
2279 covTr,
2280 (Short_t)esdTrack->GetSign(),
2281 esdTrack->GetITSClusterMap(),
2282 pid,
2283 vcascade,
2284 kTRUE, // check if this is right
2285 kFALSE, // check if this is right
2286 AliAODTrack::kSecondary)
2287 );
2288 aodTrack->ConvertAliPIDtoAODPID();
2289 }
2290 else {
2291 cerr << "Error: event " << iEvent << " cascade " << nCascade
2292 << " track " << bachelor << " has already been used!" << endl;
2293 }
2294
2295 // Add the primary track of the cascade (if any)
2296
2297 } // end of the loop on cascades
2298
2299 // V0s
2300
2301 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2302
2303 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2304
2305 AliESDv0 *v0 = esd->GetV0(nV0);
2306
2307 v0->GetXYZ(pos[0], pos[1], pos[2]);
2308 v0->GetPosCov(covVtx);
2309
2310 AliAODVertex * vV0 =
2311 new(vertices[jVertices++]) AliAODVertex(pos,
2312 covVtx,
2313 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2314 primary,
2315 AliAODVertex::kV0);
2316 primary->AddDaughter(vV0);
2317
2318 Int_t posFromV0 = v0->GetPindex();
2319 Int_t negFromV0 = v0->GetNindex();
2320
2321 // Add the positive tracks from the V0
2322
2323 if (!usedTrack[posFromV0]) {
2324
2325 usedTrack[posFromV0] = kTRUE;
2326
2327 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2328 esdTrack->GetPxPyPz(p);
2329 esdTrack->GetXYZ(pos);
2330 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2331 esdTrack->GetESDpid(pid);
2332
2333 vV0->AddDaughter(aodTrack =
2334 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2335 esdTrack->GetLabel(),
2336 p,
2337 kTRUE,
2338 pos,
2339 kFALSE,
2340 covTr,
2341 (Short_t)esdTrack->GetSign(),
2342 esdTrack->GetITSClusterMap(),
2343 pid,
2344 vV0,
2345 kTRUE, // check if this is right
2346 kFALSE, // check if this is right
2347 AliAODTrack::kSecondary)
2348 );
2349 aodTrack->ConvertAliPIDtoAODPID();
2350 }
2351 else {
2352 cerr << "Error: event " << iEvent << " V0 " << nV0
2353 << " track " << posFromV0 << " has already been used!" << endl;
2354 }
a7807689 2355
f29f1726 2356 // Add the negative tracks from the V0
2357
2358 if (!usedTrack[negFromV0]) {
2359
2360 usedTrack[negFromV0] = kTRUE;
2361
2362 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2363 esdTrack->GetPxPyPz(p);
2364 esdTrack->GetXYZ(pos);
2365 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2366 esdTrack->GetESDpid(pid);
2367
2368 vV0->AddDaughter(aodTrack =
2369 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2370 esdTrack->GetLabel(),
2371 p,
2372 kTRUE,
2373 pos,
2374 kFALSE,
2375 covTr,
2376 (Short_t)esdTrack->GetSign(),
2377 esdTrack->GetITSClusterMap(),
2378 pid,
2379 vV0,
2380 kTRUE, // check if this is right
2381 kFALSE, // check if this is right
2382 AliAODTrack::kSecondary)
2383 );
2384 aodTrack->ConvertAliPIDtoAODPID();
2385 }
2386 else {
2387 cerr << "Error: event " << iEvent << " V0 " << nV0
2388 << " track " << negFromV0 << " has already been used!" << endl;
2389 }
2390
2391 } // end of the loop on V0s
2392
2393 // Kinks: it is a big mess the access to the information in the kinks
2394 // The loop is on the tracks in order to find the mother and daugther of each kink
2395
2396
2397 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2398
2399
2400 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2401
2402 Int_t ikink = esdTrack->GetKinkIndex(0);
2403
2404 if (ikink) {
2405 // Negative kink index: mother, positive: daughter
2406
2407 // Search for the second track of the kink
2408
2409 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2410
2411 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2412
2413 Int_t jkink = esdTrack1->GetKinkIndex(0);
2414
2415 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2416
2417 // The two tracks are from the same kink
2418
2419 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2420
2421 Int_t imother = -1;
2422 Int_t idaughter = -1;
2423
2424 if (ikink<0 && jkink>0) {
2425
2426 imother = iTrack;
2427 idaughter = jTrack;
2428 }
2429 else if (ikink>0 && jkink<0) {
2430
2431 imother = jTrack;
2432 idaughter = iTrack;
2433 }
2434 else {
2435 cerr << "Error: Wrong combination of kink indexes: "
2436 << ikink << " " << jkink << endl;
2437 continue;
2438 }
2439
2440 // Add the mother track
2441
2442 AliAODTrack * mother = NULL;
2443
2444 if (!usedTrack[imother]) {
2445
2446 usedTrack[imother] = kTRUE;
2447
2448 AliESDtrack *esdTrack = esd->GetTrack(imother);
2449 esdTrack->GetPxPyPz(p);
2450 esdTrack->GetXYZ(pos);
2451 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2452 esdTrack->GetESDpid(pid);
2453
2454 mother =
2455 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2456 esdTrack->GetLabel(),
2457 p,
2458 kTRUE,
2459 pos,
2460 kFALSE,
2461 covTr,
2462 (Short_t)esdTrack->GetSign(),
2463 esdTrack->GetITSClusterMap(),
2464 pid,
2465 primary,
2466 kTRUE, // check if this is right
2467 kTRUE, // check if this is right
2468 AliAODTrack::kPrimary);
2469 primary->AddDaughter(mother);
2470 mother->ConvertAliPIDtoAODPID();
2471 }
2472 else {
2473 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2474 << " track " << imother << " has already been used!" << endl;
2475 }
2476
2477 // Add the kink vertex
2478 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2479
2480 AliAODVertex * vkink =
2481 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2482 NULL,
2483 0.,
2484 mother,
2485 AliAODVertex::kKink);
2486 // Add the daughter track
2487
2488 AliAODTrack * daughter = NULL;
2489
2490 if (!usedTrack[idaughter]) {
2491
2492 usedTrack[idaughter] = kTRUE;
2493
2494 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2495 esdTrack->GetPxPyPz(p);
2496 esdTrack->GetXYZ(pos);
2497 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2498 esdTrack->GetESDpid(pid);
2499
2500 daughter =
2501 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2502 esdTrack->GetLabel(),
2503 p,
2504 kTRUE,
2505 pos,
2506 kFALSE,
2507 covTr,
2508 (Short_t)esdTrack->GetSign(),
2509 esdTrack->GetITSClusterMap(),
2510 pid,
2511 vkink,
2512 kTRUE, // check if this is right
2513 kTRUE, // check if this is right
2514 AliAODTrack::kPrimary);
2515 vkink->AddDaughter(daughter);
2516 daughter->ConvertAliPIDtoAODPID();
2517 }
2518 else {
2519 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2520 << " track " << idaughter << " has already been used!" << endl;
2521 }
2522
2523
2524 }
2525 }
2526
2527 }
2528
2529 }
2530
2531
2532 // Tracks (primary and orphan)
2533
2534 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2535
2536
2537 if (usedTrack[nTrack]) continue;
2538
2539 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2540 esdTrack->GetPxPyPz(p);
2541 esdTrack->GetXYZ(pos);
2542 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2543 esdTrack->GetESDpid(pid);
2544
2545 Float_t impactXY, impactZ;
2546
2547 esdTrack->GetImpactParameters(impactXY,impactZ);
2548
2549 if (impactXY<3) {
2550 // track inside the beam pipe
2551
2552 primary->AddDaughter(aodTrack =
2553 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2554 esdTrack->GetLabel(),
2555 p,
2556 kTRUE,
2557 pos,
2558 kFALSE,
2559 covTr,
2560 (Short_t)esdTrack->GetSign(),
2561 esdTrack->GetITSClusterMap(),
2562 pid,
2563 primary,
2564 kTRUE, // check if this is right
2565 kTRUE, // check if this is right
2566 AliAODTrack::kPrimary)
2567 );
2568 aodTrack->ConvertAliPIDtoAODPID();
2569 }
2570 else {
2571 // outside the beam pipe: orphan track
2572 aodTrack =
2573 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2574 esdTrack->GetLabel(),
2575 p,
2576 kTRUE,
2577 pos,
2578 kFALSE,
2579 covTr,
2580 (Short_t)esdTrack->GetSign(),
2581 esdTrack->GetITSClusterMap(),
2582 pid,
2583 NULL,
2584 kFALSE, // check if this is right
2585 kFALSE, // check if this is right
2586 AliAODTrack::kOrphan);
2587 aodTrack->ConvertAliPIDtoAODPID();
2588 }
2589 } // end of loop on tracks
2590
2591 // muon tracks
2592 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2593 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2594
2595 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2596 p[0] = esdMuTrack->Px();
2597 p[1] = esdMuTrack->Py();
2598 p[2] = esdMuTrack->Pz();
2599 pos[0] = primary->GetX();
2600 pos[1] = primary->GetY();
2601 pos[2] = primary->GetZ();
2602
2603 // has to be changed once the muon pid is provided by the ESD
2604 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2605
2606 primary->AddDaughter(
2607 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2608 0, // no label provided
2609 p,
2610 kTRUE,
2611 pos,
2612 kFALSE,
2613 NULL, // no covariance matrix provided
2614 (Short_t)-99, // no charge provided
2615 0, // no ITSClusterMap
2616 pid,
2617 primary,
2618 kTRUE, // check if this is right
2619 kTRUE, // not used for vertex fit
2620 AliAODTrack::kPrimary)
2621 );
2622 }
2623
2624 // Access to the AOD container of clusters
2625 TClonesArray &clusters = *(aod->GetClusters());
2626 Int_t jClusters=0;
2627
2628 // Calo Clusters
2629 Int_t nClusters = esd->GetNumberOfCaloClusters();
2630
2631 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2632
2633 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2634
2635 Int_t id = cluster->GetID();
2636 Int_t label = -1;
53ec9628 2637 Float_t energy = cluster->E();
2638 cluster->GetPosition(posF);
f29f1726 2639 AliAODVertex *prodVertex = primary;
2640 AliAODTrack *primTrack = NULL;
2641 Char_t ttype=AliAODCluster::kUndef;
2642
2643 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2644 else if (cluster->IsEMCAL()) {
2645
2646 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2647 ttype = AliAODCluster::kEMCALPseudoCluster;
2648 else
2649 ttype = AliAODCluster::kEMCALClusterv1;
2650
2651 }
2652
2653 new(clusters[jClusters++]) AliAODCluster(id,
2654 label,
2655 energy,
2656 pos,
2657 NULL, // no covariance matrix provided
2658 NULL, // no pid for clusters provided
2659 prodVertex,
2660 primTrack,
2661 ttype);
2662
2663 } // end of loop on calo clusters
2664
2665 delete [] usedTrack;
2666 delete [] usedV0;
2667 delete [] usedKink;
2668
2669 // fill the tree for this event
2670 aodTree->Fill();
2671 } // end of event loop
2672
2673 aodTree->GetUserInfo()->Add(aod);
2674
2675 // close ESD file
2676 esdFile->Close();
2677
2678 // write the tree to the specified file
2679 aodFile = aodTree->GetCurrentFile();
2680 aodFile->cd();
2681 aodTree->Write();
2682
a7807689 2683 return;
2684}
2685
98937d93 2686void AliReconstruction::WriteAlignmentData(AliESD* esd)
2687{
2688 // Write space-points which are then used in the alignment procedures
2689 // For the moment only ITS, TRD and TPC
2690
2691 // Load TOF clusters
d528ee75 2692 if (fTracker[3]){
2693 fLoader[3]->LoadRecPoints("read");
2694 TTree* tree = fLoader[3]->TreeR();
2695 if (!tree) {
2696 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2697 return;
2698 }
2699 fTracker[3]->LoadClusters(tree);
98937d93 2700 }
98937d93 2701 Int_t ntracks = esd->GetNumberOfTracks();
2702 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2703 {
2704 AliESDtrack *track = esd->GetTrack(itrack);
2705 Int_t nsp = 0;
ef7253ac 2706 Int_t idx[200];
98937d93 2707 for (Int_t iDet = 3; iDet >= 0; iDet--)
2708 nsp += track->GetNcls(iDet);
2709 if (nsp) {
2710 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2711 track->SetTrackPointArray(sp);
2712 Int_t isptrack = 0;
2713 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2714 AliTracker *tracker = fTracker[iDet];
2715 if (!tracker) continue;
2716 Int_t nspdet = track->GetNcls(iDet);
98937d93 2717 if (nspdet <= 0) continue;
2718 track->GetClusters(iDet,idx);
2719 AliTrackPoint p;
2720 Int_t isp = 0;
2721 Int_t isp2 = 0;
2722 while (isp < nspdet) {
2723 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 2724 const Int_t kNTPCmax = 159;
2725 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2726 if (!isvalid) continue;
2727 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2728 }
98937d93 2729 }
2730 }
2731 }
d528ee75 2732 if (fTracker[3]){
2733 fTracker[3]->UnloadClusters();
2734 fLoader[3]->UnloadRecPoints();
2735 }
98937d93 2736}
2e3550da 2737
2738//_____________________________________________________________________________
2739void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2740{
2741 // The method reads the raw-data error log
2742 // accumulated within the rawReader.
2743 // It extracts the raw-data errors related to
2744 // the current event and stores them into
2745 // a TClonesArray inside the esd object.
2746
2747 if (!fRawReader) return;
2748
2749 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2750
2751 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2752 if (!log) continue;
2753 if (iEvent != log->GetEventNumber()) continue;
2754
2755 esd->AddRawDataErrorLog(log);
2756 }
2757
2758}
46698ae4 2759
2760TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
b545009a 2761 // Dump a file content into a char in TNamed
46698ae4 2762 ifstream in;
2763 in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2764 Int_t kBytes = (Int_t)in.tellg();
2765 printf("Size: %d \n",kBytes);
2766 TNamed *fn = 0;
2767 if(in.good()){
2768 char* memblock = new char [kBytes];
2769 in.seekg (0, ios::beg);
2770 in.read (memblock, kBytes);
2771 in.close();
2772 TString fData(memblock,kBytes);
2773 fn = new TNamed(fName,fData);
2774 printf("fData Size: %d \n",fData.Sizeof());
2775 printf("fName Size: %d \n",fName.Sizeof());
2776 printf("fn Size: %d \n",fn->Sizeof());
2777 delete[] memblock;
2778 }
2779 else{
2780 AliInfo(Form("Could not Open %s\n",fPath.Data()));
2781 }
2782
2783 return fn;
2784}
2785
2786void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
46698ae4 2787 // This is not really needed in AliReconstruction at the moment
2788 // but can serve as a template
2789
2790 TList *fList = fTree->GetUserInfo();
2791 TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2792 printf("fn Size: %d \n",fn->Sizeof());
2793
2794 TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2795 const char* cdata = fn->GetTitle();
2796 printf("fTmp Size %d\n",fTmp.Sizeof());
2797
2798 int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2799 printf("calculated size %d\n",size);
2800 ofstream out(fName.Data(),ios::out | ios::binary);
2801 out.write(cdata,size);
2802 out.close();
2803
2804}
2805