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