]> 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
Reconstruction QA by Sylwester
[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 {
944 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 945 vertex->SetName("reconstructed");
2257f27e 946 }
947
948 } else {
815c2b38 949 AliInfo("getting the primary vertex from MC");
2257f27e 950 vertex = new AliESDVertex(vtxPos, vtxErr);
951 }
952
953 if (vertex) {
954 vertex->GetXYZ(vtxPos);
955 vertex->GetSigmaXYZ(vtxErr);
956 } else {
815c2b38 957 AliWarning("no vertex reconstructed");
2257f27e 958 vertex = new AliESDVertex(vtxPos, vtxErr);
959 }
960 esd->SetVertex(vertex);
32e449be 961 // if SPD multiplicity has been determined, it is stored in the ESD
25be1e5c 962 AliMultiplicity *mult = fVertexer->GetMultiplicity();
32e449be 963 if(mult)esd->SetMultiplicity(mult);
964
b8cd5251 965 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
966 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
967 }
2257f27e 968 delete vertex;
969
5f8272e1 970 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
971 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 972
973 return kTRUE;
974}
975
1f46a9ae 976//_____________________________________________________________________________
977Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
978{
979// run the HLT barrel tracking
980
981 TStopwatch stopwatch;
982 stopwatch.Start();
983
984 if (!fRunLoader) {
985 AliError("Missing runLoader!");
986 return kFALSE;
987 }
988
989 AliInfo("running HLT tracking");
990
991 // Get a pointer to the HLT reconstructor
992 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
993 if (!reconstructor) return kFALSE;
994
995 // TPC + ITS
996 for (Int_t iDet = 1; iDet >= 0; iDet--) {
997 TString detName = fgkDetectorName[iDet];
998 AliDebug(1, Form("%s HLT tracking", detName.Data()));
999 reconstructor->SetOption(detName.Data());
1000 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1001 if (!tracker) {
1002 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1003 if (fStopOnError) return kFALSE;
9dcc06e1 1004 continue;
1f46a9ae 1005 }
1006 Double_t vtxPos[3];
1007 Double_t vtxErr[3]={0.005,0.005,0.010};
1008 const AliESDVertex *vertex = esd->GetVertex();
1009 vertex->GetXYZ(vtxPos);
1010 tracker->SetVertex(vtxPos,vtxErr);
1011 if(iDet != 1) {
1012 fLoader[iDet]->LoadRecPoints("read");
1013 TTree* tree = fLoader[iDet]->TreeR();
1014 if (!tree) {
1015 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1016 return kFALSE;
1017 }
1018 tracker->LoadClusters(tree);
1019 }
1020 if (tracker->Clusters2Tracks(esd) != 0) {
1021 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1022 return kFALSE;
1023 }
1024 if(iDet != 1) {
1025 tracker->UnloadClusters();
1026 }
1027 delete tracker;
1028 }
1029
5f8272e1 1030 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1031 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 1032
1033 return kTRUE;
1034}
1035
e66fbafb 1036//_____________________________________________________________________________
1037Bool_t AliReconstruction::RunMuonTracking()
1038{
1039// run the muon spectrometer tracking
1040
1041 TStopwatch stopwatch;
1042 stopwatch.Start();
1043
1044 if (!fRunLoader) {
1045 AliError("Missing runLoader!");
1046 return kFALSE;
1047 }
1048 Int_t iDet = 7; // for MUON
1049
1050 AliInfo("is running...");
1051
1052 // Get a pointer to the MUON reconstructor
1053 AliReconstructor *reconstructor = GetReconstructor(iDet);
1054 if (!reconstructor) return kFALSE;
1055
1056
1057 TString detName = fgkDetectorName[iDet];
1058 AliDebug(1, Form("%s tracking", detName.Data()));
1059 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1060 if (!tracker) {
1061 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1062 return kFALSE;
1063 }
1064
1065 // create Tracks
1066 fLoader[iDet]->LoadTracks("update");
1067 fLoader[iDet]->CleanTracks();
1068 fLoader[iDet]->MakeTracksContainer();
1069
1070 // read RecPoints
1071 fLoader[iDet]->LoadRecPoints("read");
1072
1073 if (!tracker->Clusters2Tracks(0x0)) {
1074 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1075 return kFALSE;
1076 }
1077 fLoader[iDet]->UnloadRecPoints();
1078
1079 fLoader[iDet]->WriteTracks("OVERWRITE");
1080 fLoader[iDet]->UnloadTracks();
1081
1082 delete tracker;
1083
1084
1085 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1086 stopwatch.RealTime(),stopwatch.CpuTime()));
1087
1088 return kTRUE;
1089}
1090
1091
2257f27e 1092//_____________________________________________________________________________
1093Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1094{
1095// run the barrel tracking
1096
1097 TStopwatch stopwatch;
1098 stopwatch.Start();
24f7a148 1099
815c2b38 1100 AliInfo("running tracking");
596a855f 1101
91b876d1 1102 //Fill the ESD with the T0 info (will be used by the TOF)
1f4331b3 1103 if (fReconstructor[11])
1104 GetReconstructor(11)->FillESD(fRunLoader, esd);
91b876d1 1105
b8cd5251 1106 // pass 1: TPC + ITS inwards
1107 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1108 if (!fTracker[iDet]) continue;
1109 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 1110
b8cd5251 1111 // load clusters
1112 fLoader[iDet]->LoadRecPoints("read");
1113 TTree* tree = fLoader[iDet]->TreeR();
1114 if (!tree) {
1115 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1116 return kFALSE;
1117 }
b8cd5251 1118 fTracker[iDet]->LoadClusters(tree);
1119
1120 // run tracking
1121 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1122 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 1123 return kFALSE;
1124 }
b8cd5251 1125 if (fCheckPointLevel > 1) {
1126 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1127 }
878e1fe1 1128 // preliminary PID in TPC needed by the ITS tracker
1129 if (iDet == 1) {
1130 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 1131 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 1132 AliESDpid::MakePID(esd);
1133 }
b8cd5251 1134 }
596a855f 1135
b8cd5251 1136 // pass 2: ALL backwards
1137 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1138 if (!fTracker[iDet]) continue;
1139 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1140
1141 // load clusters
1142 if (iDet > 1) { // all except ITS, TPC
1143 TTree* tree = NULL;
7b61cd9c 1144 fLoader[iDet]->LoadRecPoints("read");
1145 tree = fLoader[iDet]->TreeR();
b8cd5251 1146 if (!tree) {
1147 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1148 return kFALSE;
1149 }
b8cd5251 1150 fTracker[iDet]->LoadClusters(tree);
1151 }
24f7a148 1152
b8cd5251 1153 // run tracking
1154 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1155 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1156 return kFALSE;
1157 }
1158 if (fCheckPointLevel > 1) {
1159 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1160 }
24f7a148 1161
b8cd5251 1162 // unload clusters
1163 if (iDet > 2) { // all except ITS, TPC, TRD
1164 fTracker[iDet]->UnloadClusters();
7b61cd9c 1165 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1166 }
8f37df88 1167 // updated PID in TPC needed by the ITS tracker -MI
1168 if (iDet == 1) {
1169 GetReconstructor(1)->FillESD(fRunLoader, esd);
1170 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1171 AliESDpid::MakePID(esd);
1172 }
b8cd5251 1173 }
596a855f 1174
98937d93 1175 // write space-points to the ESD in case alignment data output
1176 // is switched on
1177 if (fWriteAlignmentData)
1178 WriteAlignmentData(esd);
1179
b8cd5251 1180 // pass 3: TRD + TPC + ITS refit inwards
1181 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1182 if (!fTracker[iDet]) continue;
1183 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 1184
b8cd5251 1185 // run tracking
1186 if (fTracker[iDet]->RefitInward(esd) != 0) {
1187 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1188 return kFALSE;
1189 }
1190 if (fCheckPointLevel > 1) {
1191 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1192 }
596a855f 1193
b8cd5251 1194 // unload clusters
1195 fTracker[iDet]->UnloadClusters();
1196 fLoader[iDet]->UnloadRecPoints();
1197 }
ff8bb5ae 1198 //
1199 // Propagate track to the vertex - if not done by ITS
1200 //
1201 Int_t ntracks = esd->GetNumberOfTracks();
1202 for (Int_t itrack=0; itrack<ntracks; itrack++){
1203 const Double_t kRadius = 3; // beam pipe radius
1204 const Double_t kMaxStep = 5; // max step
1205 const Double_t kMaxD = 123456; // max distance to prim vertex
1206 Double_t fieldZ = AliTracker::GetBz(); //
1207 AliESDtrack * track = esd->GetTrack(itrack);
1208 if (!track) continue;
1209 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
a28195f7 1210 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
ff8bb5ae 1211 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1212 }
1213
5f8272e1 1214 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1215 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1216
596a855f 1217 return kTRUE;
1218}
1219
1220//_____________________________________________________________________________
24f7a148 1221Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1222{
1223// fill the event summary data
1224
030b532d 1225 TStopwatch stopwatch;
1226 stopwatch.Start();
815c2b38 1227 AliInfo("filling ESD");
030b532d 1228
596a855f 1229 TString detStr = detectors;
b8cd5251 1230 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1231 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1232 AliReconstructor* reconstructor = GetReconstructor(iDet);
1233 if (!reconstructor) continue;
1234
1235 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1236 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 1237 TTree* clustersTree = NULL;
1238 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1239 fLoader[iDet]->LoadRecPoints("read");
1240 clustersTree = fLoader[iDet]->TreeR();
1241 if (!clustersTree) {
1242 AliError(Form("Can't get the %s clusters tree",
1243 fgkDetectorName[iDet]));
1244 if (fStopOnError) return kFALSE;
1245 }
1246 }
1247 if (fRawReader && !reconstructor->HasDigitConversion()) {
1248 reconstructor->FillESD(fRawReader, clustersTree, esd);
1249 } else {
1250 TTree* digitsTree = NULL;
1251 if (fLoader[iDet]) {
1252 fLoader[iDet]->LoadDigits("read");
1253 digitsTree = fLoader[iDet]->TreeD();
1254 if (!digitsTree) {
1255 AliError(Form("Can't get the %s digits tree",
1256 fgkDetectorName[iDet]));
1257 if (fStopOnError) return kFALSE;
1258 }
1259 }
1260 reconstructor->FillESD(digitsTree, clustersTree, esd);
1261 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1262 }
1263 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1264 fLoader[iDet]->UnloadRecPoints();
1265 }
1266
b8cd5251 1267 if (fRawReader) {
1268 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1269 } else {
1270 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1271 }
b8cd5251 1272 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1273 }
1274 }
1275
1276 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1277 AliError(Form("the following detectors were not found: %s",
1278 detStr.Data()));
596a855f 1279 if (fStopOnError) return kFALSE;
1280 }
1281
5f8272e1 1282 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1283 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1284
596a855f 1285 return kTRUE;
1286}
1287
b647652d 1288//_____________________________________________________________________________
1289Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1290{
1291 // Reads the trigger decision which is
1292 // stored in Trigger.root file and fills
1293 // the corresponding esd entries
1294
1295 AliInfo("Filling trigger information into the ESD");
1296
1297 if (fRawReader) {
1298 AliCTPRawStream input(fRawReader);
1299 if (!input.Next()) {
1300 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1301 return kFALSE;
1302 }
1303 esd->SetTriggerMask(input.GetClassMask());
1304 esd->SetTriggerCluster(input.GetClusterMask());
1305 }
1306 else {
1307 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1308 if (runloader) {
1309 if (!runloader->LoadTrigger()) {
1310 AliCentralTrigger *aCTP = runloader->GetTrigger();
1311 esd->SetTriggerMask(aCTP->GetClassMask());
1312 esd->SetTriggerCluster(aCTP->GetClusterMask());
1313 }
1314 else {
1315 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1316 return kFALSE;
1317 }
1318 }
1319 else {
1320 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1321 return kFALSE;
1322 }
1323 }
1324
1325 return kTRUE;
1326}
596a855f 1327
001397cd 1328
1329
1330
1331
1332//_____________________________________________________________________________
1333Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1334{
1335 //
1336 // Filling information from RawReader Header
1337 //
1338
1339 AliInfo("Filling information from RawReader Header");
31fd97b2 1340 esd->SetBunchCrossNumber(0);
1341 esd->SetOrbitNumber(0);
9bcc1e45 1342 esd->SetPeriodNumber(0);
001397cd 1343 esd->SetTimeStamp(0);
1344 esd->SetEventType(0);
1345 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1346 if (eventHeader){
9bcc1e45 1347
1348 const UInt_t *id = eventHeader->GetP("Id");
1349 esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1350 esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1351 esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1352
001397cd 1353 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
31fd97b2 1354 esd->SetEventType((eventHeader->Get("Type")));
001397cd 1355 }
1356
1357 return kTRUE;
1358}
1359
1360
596a855f 1361//_____________________________________________________________________________
1362Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1363{
1364// check whether detName is contained in detectors
1365// if yes, it is removed from detectors
1366
1367 // check if all detectors are selected
1368 if ((detectors.CompareTo("ALL") == 0) ||
1369 detectors.BeginsWith("ALL ") ||
1370 detectors.EndsWith(" ALL") ||
1371 detectors.Contains(" ALL ")) {
1372 detectors = "ALL";
1373 return kTRUE;
1374 }
1375
1376 // search for the given detector
1377 Bool_t result = kFALSE;
1378 if ((detectors.CompareTo(detName) == 0) ||
1379 detectors.BeginsWith(detName+" ") ||
1380 detectors.EndsWith(" "+detName) ||
1381 detectors.Contains(" "+detName+" ")) {
1382 detectors.ReplaceAll(detName, "");
1383 result = kTRUE;
1384 }
1385
1386 // clean up the detectors string
1387 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1388 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1389 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1390
1391 return result;
1392}
e583c30d 1393
f08fc9f5 1394//_____________________________________________________________________________
1395Bool_t AliReconstruction::InitRunLoader()
1396{
1397// get or create the run loader
1398
1399 if (gAlice) delete gAlice;
1400 gAlice = NULL;
1401
b26c3770 1402 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1403 // load all base libraries to get the loader classes
1404 TString libs = gSystem->GetLibraries();
1405 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1406 TString detName = fgkDetectorName[iDet];
1407 if (detName == "HLT") continue;
1408 if (libs.Contains("lib" + detName + "base.so")) continue;
1409 gSystem->Load("lib" + detName + "base.so");
1410 }
f08fc9f5 1411 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1412 if (!fRunLoader) {
1413 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1414 CleanUp();
1415 return kFALSE;
1416 }
b26c3770 1417 fRunLoader->CdGAFile();
1418 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1419 if (fRunLoader->LoadgAlice() == 0) {
1420 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1421 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1422 }
f08fc9f5 1423 }
1424 if (!gAlice && !fRawReader) {
1425 AliError(Form("no gAlice object found in file %s",
1426 fGAliceFileName.Data()));
1427 CleanUp();
1428 return kFALSE;
1429 }
1430
1431 } else { // galice.root does not exist
1432 if (!fRawReader) {
1433 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1434 CleanUp();
1435 return kFALSE;
1436 }
1437 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1438 AliConfig::GetDefaultEventFolderName(),
1439 "recreate");
1440 if (!fRunLoader) {
1441 AliError(Form("could not create run loader in file %s",
1442 fGAliceFileName.Data()));
1443 CleanUp();
1444 return kFALSE;
1445 }
1446 fRunLoader->MakeTree("E");
1447 Int_t iEvent = 0;
1448 while (fRawReader->NextEvent()) {
1449 fRunLoader->SetEventNumber(iEvent);
1450 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1451 iEvent, iEvent);
1452 fRunLoader->MakeTree("H");
1453 fRunLoader->TreeE()->Fill();
1454 iEvent++;
1455 }
1456 fRawReader->RewindEvents();
973388c2 1457 if (fNumberOfEventsPerFile > 0)
1458 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1459 else
1460 fRunLoader->SetNumberOfEventsPerFile(iEvent);
f08fc9f5 1461 fRunLoader->WriteHeader("OVERWRITE");
1462 fRunLoader->CdGAFile();
1463 fRunLoader->Write(0, TObject::kOverwrite);
1464// AliTracker::SetFieldMap(???);
1465 }
1466
1467 return kTRUE;
1468}
1469
c757bafd 1470//_____________________________________________________________________________
b8cd5251 1471AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1472{
f08fc9f5 1473// get the reconstructor object and the loader for a detector
c757bafd 1474
b8cd5251 1475 if (fReconstructor[iDet]) return fReconstructor[iDet];
1476
1477 // load the reconstructor object
1478 TPluginManager* pluginManager = gROOT->GetPluginManager();
1479 TString detName = fgkDetectorName[iDet];
1480 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1481 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1482
1483 if (detName == "HLT") {
1484 if (!gROOT->GetClass("AliLevel3")) {
4aa41877 1485 gSystem->Load("libAliHLTSrc.so");
1486 gSystem->Load("libAliHLTMisc.so");
1487 gSystem->Load("libAliHLTHough.so");
1488 gSystem->Load("libAliHLTComp.so");
b8cd5251 1489 }
1490 }
1491
1492 AliReconstructor* reconstructor = NULL;
1493 // first check if a plugin is defined for the reconstructor
1494 TPluginHandler* pluginHandler =
1495 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1496 // if not, add a plugin for it
1497 if (!pluginHandler) {
b8cd5251 1498 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1499 TString libs = gSystem->GetLibraries();
1500 if (libs.Contains("lib" + detName + "base.so") ||
1501 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1502 pluginManager->AddHandler("AliReconstructor", detName,
1503 recName, detName + "rec", recName + "()");
1504 } else {
1505 pluginManager->AddHandler("AliReconstructor", detName,
1506 recName, detName, recName + "()");
c757bafd 1507 }
b8cd5251 1508 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1509 }
1510 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1511 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1512 }
b8cd5251 1513 if (reconstructor) {
1514 TObject* obj = fOptions.FindObject(detName.Data());
1515 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1516 reconstructor->Init(fRunLoader);
b8cd5251 1517 fReconstructor[iDet] = reconstructor;
1518 }
1519
f08fc9f5 1520 // get or create the loader
1521 if (detName != "HLT") {
1522 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1523 if (!fLoader[iDet]) {
1524 AliConfig::Instance()
1525 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1526 detName, detName);
1527 // first check if a plugin is defined for the loader
bb0901a4 1528 pluginHandler =
f08fc9f5 1529 pluginManager->FindHandler("AliLoader", detName);
1530 // if not, add a plugin for it
1531 if (!pluginHandler) {
1532 TString loaderName = "Ali" + detName + "Loader";
1533 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1534 pluginManager->AddHandler("AliLoader", detName,
1535 loaderName, detName + "base",
1536 loaderName + "(const char*, TFolder*)");
1537 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1538 }
1539 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1540 fLoader[iDet] =
1541 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1542 fRunLoader->GetEventFolder());
1543 }
1544 if (!fLoader[iDet]) { // use default loader
1545 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1546 }
1547 if (!fLoader[iDet]) {
1548 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1549 if (fStopOnError) return NULL;
f08fc9f5 1550 } else {
1551 fRunLoader->AddLoader(fLoader[iDet]);
1552 fRunLoader->CdGAFile();
1553 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1554 fRunLoader->Write(0, TObject::kOverwrite);
1555 }
1556 }
1557 }
1558
b8cd5251 1559 return reconstructor;
c757bafd 1560}
1561
2257f27e 1562//_____________________________________________________________________________
1563Bool_t AliReconstruction::CreateVertexer()
1564{
1565// create the vertexer
1566
b8cd5251 1567 fVertexer = NULL;
1568 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1569 if (itsReconstructor) {
b8cd5251 1570 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1571 }
b8cd5251 1572 if (!fVertexer) {
815c2b38 1573 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1574 if (fStopOnError) return kFALSE;
1575 }
1576
1577 return kTRUE;
1578}
1579
24f7a148 1580//_____________________________________________________________________________
b8cd5251 1581Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1582{
f08fc9f5 1583// create the trackers
24f7a148 1584
b8cd5251 1585 TString detStr = detectors;
1586 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1587 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1588 AliReconstructor* reconstructor = GetReconstructor(iDet);
1589 if (!reconstructor) continue;
1590 TString detName = fgkDetectorName[iDet];
1f46a9ae 1591 if (detName == "HLT") {
1592 fRunHLTTracking = kTRUE;
1593 continue;
1594 }
e66fbafb 1595 if (detName == "MUON") {
1596 fRunMuonTracking = kTRUE;
1597 continue;
1598 }
1599
f08fc9f5 1600
1601 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1602 if (!fTracker[iDet] && (iDet < 7)) {
1603 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1604 if (fStopOnError) return kFALSE;
1605 }
1606 }
1607
24f7a148 1608 return kTRUE;
1609}
1610
e583c30d 1611//_____________________________________________________________________________
b26c3770 1612void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1613{
1614// delete trackers and the run loader and close and delete the file
1615
b8cd5251 1616 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1617 delete fReconstructor[iDet];
1618 fReconstructor[iDet] = NULL;
1619 fLoader[iDet] = NULL;
1620 delete fTracker[iDet];
1621 fTracker[iDet] = NULL;
1622 }
1623 delete fVertexer;
1624 fVertexer = NULL;
9178838a 1625 delete fDiamondProfile;
1626 fDiamondProfile = NULL;
e583c30d 1627
1628 delete fRunLoader;
1629 fRunLoader = NULL;
b649205a 1630 delete fRawReader;
1631 fRawReader = NULL;
e583c30d 1632
1633 if (file) {
1634 file->Close();
1635 delete file;
1636 }
b26c3770 1637
1638 if (fileOld) {
1639 fileOld->Close();
1640 delete fileOld;
1641 gSystem->Unlink("AliESDs.old.root");
1642 }
e583c30d 1643}
24f7a148 1644
1645
1646//_____________________________________________________________________________
1647Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1648{
1649// read the ESD event from a file
1650
1651 if (!esd) return kFALSE;
1652 char fileName[256];
1653 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1654 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1655 if (gSystem->AccessPathName(fileName)) return kFALSE;
1656
f3a97c86 1657 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1658 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1659 TFile* file = TFile::Open(fileName);
1660 if (!file || !file->IsOpen()) {
815c2b38 1661 AliError(Form("opening %s failed", fileName));
24f7a148 1662 delete file;
1663 return kFALSE;
1664 }
1665
1666 gROOT->cd();
1667 delete esd;
1668 esd = (AliESD*) file->Get("ESD");
1669 file->Close();
1670 delete file;
1671 return kTRUE;
1672}
1673
1674//_____________________________________________________________________________
1675void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1676{
1677// write the ESD event to a file
1678
1679 if (!esd) return;
1680 char fileName[256];
1681 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1682 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1683
815c2b38 1684 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1685 TFile* file = TFile::Open(fileName, "recreate");
1686 if (!file || !file->IsOpen()) {
815c2b38 1687 AliError(Form("opening %s failed", fileName));
24f7a148 1688 } else {
1689 esd->Write("ESD");
1690 file->Close();
1691 }
1692 delete file;
1693}
f3a97c86 1694
1695
1696
1697
1698//_____________________________________________________________________________
1699void AliReconstruction::CreateTag(TFile* file)
1700{
1ac7d449 1701 //GRP
1702 Float_t lhcLuminosity = 0.0;
1703 TString lhcState = "test";
1704 UInt_t detectorMask = 0;
1705
2bdb9d38 1706 /////////////
1707 //muon code//
1708 ////////////
56982dd1 1709 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1710 //Variables
56982dd1 1711 Double_t fX,fY,fZ ;
1712 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1713 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1714 Int_t fCharge;
1715 TLorentzVector fEPvector;
1716
1717 Float_t fZVertexCut = 10.0;
1718 Float_t fRhoVertexCut = 2.0;
1719
1720 Float_t fLowPtCut = 1.0;
1721 Float_t fHighPtCut = 3.0;
1722 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1723 ////////////
1724
1725 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1726
089bf903 1727 // Creates the tags for all the events in a given ESD file
4302e20f 1728 Int_t ntrack;
089bf903 1729 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1730 Int_t nPos, nNeg, nNeutr;
1731 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1732 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1733 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1734 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1735 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1736 Int_t fVertexflag;
1387d165 1737 Int_t iRunNumber = 0;
547d75a6 1738 TString fVertexName("default");
4302e20f 1739
f3a97c86 1740 AliRunTag *tag = new AliRunTag();
f3a97c86 1741 AliEventTag *evTag = new AliEventTag();
1742 TTree ttag("T","A Tree with event tags");
2bdb9d38 1743 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1744 btag->SetCompressionLevel(9);
2bdb9d38 1745
f3a97c86 1746 AliInfo(Form("Creating the tags......."));
1747
1748 if (!file || !file->IsOpen()) {
1749 AliError(Form("opening failed"));
1750 delete file;
1751 return ;
2bdb9d38 1752 }
d71aa170 1753 Int_t lastEvent = 0;
f3a97c86 1754 TTree *t = (TTree*) file->Get("esdTree");
1755 TBranch * b = t->GetBranch("ESD");
1756 AliESD *esd = 0;
1757 b->SetAddress(&esd);
1387d165 1758
d71aa170 1759 b->GetEntry(fFirstEvent);
1387d165 1760 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1761
bb0901a4 1762 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1763 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1764 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 1765 ntrack = 0;
1766 nPos = 0;
1767 nNeg = 0;
1768 nNeutr =0;
1769 nK0s = 0;
1770 nNeutrons = 0;
1771 nPi0s = 0;
1772 nGamas = 0;
1773 nProtons = 0;
1774 nKaons = 0;
1775 nPions = 0;
1776 nMuons = 0;
1777 nElectrons = 0;
1778 nCh1GeV = 0;
1779 nCh3GeV = 0;
1780 nCh10GeV = 0;
1781 nMu1GeV = 0;
1782 nMu3GeV = 0;
1783 nMu10GeV = 0;
1784 nEl1GeV = 0;
1785 nEl3GeV = 0;
1786 nEl10GeV = 0;
1787 maxPt = .0;
1788 meanPt = .0;
1789 totalP = .0;
547d75a6 1790 fVertexflag = 0;
d1a50cb5 1791
2bdb9d38 1792 b->GetEntry(iEventNumber);
1387d165 1793 iRunNumber = esd->GetRunNumber();
1794 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1795 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1796 if (!vertexIn) AliError("ESD has not defined vertex.");
1797 if (vertexIn) fVertexName = vertexIn->GetName();
1798 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1799 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1800 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1801 UInt_t status = esdTrack->GetStatus();
f3a97c86 1802
2bdb9d38 1803 //select only tracks with ITS refit
1804 if ((status&AliESDtrack::kITSrefit)==0) continue;
1805 //select only tracks with TPC refit
1806 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1807
2bdb9d38 1808 //select only tracks with the "combined PID"
1809 if ((status&AliESDtrack::kESDpid)==0) continue;
1810 Double_t p[3];
1811 esdTrack->GetPxPyPz(p);
1812 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1813 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1814 totalP += momentum;
1815 meanPt += fPt;
1816 if(fPt > maxPt) maxPt = fPt;
4302e20f 1817
2bdb9d38 1818 if(esdTrack->GetSign() > 0) {
1819 nPos++;
56982dd1 1820 if(fPt > fLowPtCut) nCh1GeV++;
1821 if(fPt > fHighPtCut) nCh3GeV++;
1822 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1823 }
1824 if(esdTrack->GetSign() < 0) {
1825 nNeg++;
56982dd1 1826 if(fPt > fLowPtCut) nCh1GeV++;
1827 if(fPt > fHighPtCut) nCh3GeV++;
1828 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1829 }
1830 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1831
2bdb9d38 1832 //PID
1833 Double_t prob[5];
1834 esdTrack->GetESDpid(prob);
4302e20f 1835
2bdb9d38 1836 Double_t rcc = 0.0;
1837 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1838 if(rcc == 0.0) continue;
1839 //Bayes' formula
1840 Double_t w[5];
1841 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1842
2bdb9d38 1843 //protons
1844 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1845 //kaons
1846 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1847 //pions
1848 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1849 //electrons
1850 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1851 nElectrons++;
56982dd1 1852 if(fPt > fLowPtCut) nEl1GeV++;
1853 if(fPt > fHighPtCut) nEl3GeV++;
1854 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1855 }
1856 ntrack++;
1857 }//track loop
1858
1859 /////////////
1860 //muon code//
1861 ////////////
1862 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1863 // loop over all reconstructed tracks (also first track of combination)
1864 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1865 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1866 if (muonTrack == 0x0) continue;
4302e20f 1867
2bdb9d38 1868 // Coordinates at vertex
56982dd1 1869 fZ = muonTrack->GetZ();
1870 fY = muonTrack->GetBendingCoor();
1871 fX = muonTrack->GetNonBendingCoor();
4302e20f 1872
56982dd1 1873 fThetaX = muonTrack->GetThetaX();
1874 fThetaY = muonTrack->GetThetaY();
4302e20f 1875
56982dd1 1876 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1877 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1878 fPxRec = fPzRec * TMath::Tan(fThetaX);
1879 fPyRec = fPzRec * TMath::Tan(fThetaY);
1880 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1881
2bdb9d38 1882 //ChiSquare of the track if needed
56982dd1 1883 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1884 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1885 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1886
2bdb9d38 1887 // total number of muons inside a vertex cut
56982dd1 1888 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1889 nMuons++;
56982dd1 1890 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1891 nMu1GeV++;
56982dd1 1892 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1893 nMu3GeV++;
56982dd1 1894 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1895 nMu10GeV++;
1896 }
1897 }
1898 }
1899 }
1900 }//muon track loop
1901
1902 // Fill the event tags
1903 if(ntrack != 0)
1904 meanPt = meanPt/ntrack;
1905
1906 evTag->SetEventId(iEventNumber+1);
547d75a6 1907 if (vertexIn) {
1908 evTag->SetVertexX(vertexIn->GetXv());
1909 evTag->SetVertexY(vertexIn->GetYv());
1910 evTag->SetVertexZ(vertexIn->GetZv());
1911 evTag->SetVertexZError(vertexIn->GetZRes());
1912 }
d1a50cb5 1913 evTag->SetVertexFlag(fVertexflag);
1914
2bdb9d38 1915 evTag->SetT0VertexZ(esd->GetT0zVertex());
1916
8bd8ac26 1917 evTag->SetTriggerMask(esd->GetTriggerMask());
1918 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1919
32a5cab4 1920 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1921 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1922 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1923 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1924 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1925 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1926
1927
1928 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1929 evTag->SetNumOfPosTracks(nPos);
1930 evTag->SetNumOfNegTracks(nNeg);
1931 evTag->SetNumOfNeutrTracks(nNeutr);
1932
1933 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1934 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1935 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1936 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1937
1938 evTag->SetNumOfProtons(nProtons);
1939 evTag->SetNumOfKaons(nKaons);
1940 evTag->SetNumOfPions(nPions);
1941 evTag->SetNumOfMuons(nMuons);
1942 evTag->SetNumOfElectrons(nElectrons);
1943 evTag->SetNumOfPhotons(nGamas);
1944 evTag->SetNumOfPi0s(nPi0s);
1945 evTag->SetNumOfNeutrons(nNeutrons);
1946 evTag->SetNumOfKaon0s(nK0s);
1947
1948 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1949 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1950 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1951 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1952 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1953 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1954 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1955 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1956 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1957
85c60a8e 1958 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1959 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1960
1961 evTag->SetTotalMomentum(totalP);
1962 evTag->SetMeanPt(meanPt);
1963 evTag->SetMaxPt(maxPt);
1964
1ac7d449 1965 tag->SetLHCTag(lhcLuminosity,lhcState);
1966 tag->SetDetectorTag(detectorMask);
1967
1387d165 1968 tag->SetRunId(iInitRunNumber);
2bdb9d38 1969 tag->AddEventTag(*evTag);
1970 }
bb0901a4 1971 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 1972 else lastEvent = fLastEvent;
f3a97c86 1973
1974 ttag.Fill();
1975 tag->Clear();
1976
1977 char fileName[256];
1978 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 1979 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 1980 AliInfo(Form("writing tags to file %s", fileName));
1981 AliDebug(1, Form("writing tags to file %s", fileName));
1982
1983 TFile* ftag = TFile::Open(fileName, "recreate");
1984 ftag->cd();
1985 ttag.Write();
1986 ftag->Close();
1987 file->cd();
1988 delete tag;
f3a97c86 1989 delete evTag;
1990}
1991
a7807689 1992//_____________________________________________________________________________
f29f1726 1993void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
a7807689 1994{
f29f1726 1995 // write all files from the given esd file to an aod file
1996
1997 // create an AliAOD object
1998 AliAODEvent *aod = new AliAODEvent();
1999 aod->CreateStdContent();
2000
2001 // go to the file
2002 aodFile->cd();
2003
2004 // create the tree
2005 TTree *aodTree = new TTree("AOD", "AliAOD tree");
2006 aodTree->Branch(aod->GetList());
2007
2008 // connect to ESD
2009 TTree *t = (TTree*) esdFile->Get("esdTree");
2010 TBranch *b = t->GetBranch("ESD");
2011 AliESD *esd = 0;
2012 b->SetAddress(&esd);
2013
2014 Int_t nEvents = b->GetEntries();
2015
2016 // set arrays and pointers
2017 Float_t posF[3];
2018 Double_t pos[3];
2019 Double_t p[3];
2020 Double_t covVtx[6];
2021 Double_t covTr[21];
2022 Double_t pid[10];
2023
2024 // loop over events and fill them
2025 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2026 b->GetEntry(iEvent);
2027
2028 // Multiplicity information needed by the header (to be revised!)
2029 Int_t nTracks = esd->GetNumberOfTracks();
2030 Int_t nPosTracks = 0;
2031 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2032 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2033
2034 // create the header
2035 aod->AddHeader(new AliAODHeader(esd->GetRunNumber(),
2036 esd->GetBunchCrossNumber(),
2037 esd->GetOrbitNumber(),
2038 esd->GetPeriodNumber(),
2039 nTracks,
2040 nPosTracks,
2041 nTracks-nPosTracks,
2042 esd->GetMagneticField(),
2043 -999., // fill muon magnetic field
2044 -999., // centrality; to be filled, still
2045 esd->GetZDCN1Energy(),
2046 esd->GetZDCP1Energy(),
2047 esd->GetZDCN2Energy(),
2048 esd->GetZDCP2Energy(),
2049 esd->GetZDCEMEnergy(),
2050 esd->GetTriggerMask(),
2051 esd->GetTriggerCluster(),
2052 esd->GetEventType()));
2053
2054 Int_t nV0s = esd->GetNumberOfV0s();
2055 Int_t nCascades = esd->GetNumberOfCascades();
2056 Int_t nKinks = esd->GetNumberOfKinks();
2057 Int_t nVertices = nV0s + nCascades + nKinks;
2058
2059 aod->ResetStd(nTracks, nVertices);
2060 AliAODTrack *aodTrack;
2061
2062
2063 // Array to take into account the tracks already added to the AOD
2064 Bool_t * usedTrack = NULL;
2065 if (nTracks>0) {
2066 usedTrack = new Bool_t[nTracks];
2067 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2068 }
2069 // Array to take into account the V0s already added to the AOD
2070 Bool_t * usedV0 = NULL;
2071 if (nV0s>0) {
2072 usedV0 = new Bool_t[nV0s];
2073 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2074 }
2075 // Array to take into account the kinks already added to the AOD
2076 Bool_t * usedKink = NULL;
2077 if (nKinks>0) {
2078 usedKink = new Bool_t[nKinks];
2079 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2080 }
2081
2082 // Access to the AOD container of vertices
2083 TClonesArray &vertices = *(aod->GetVertices());
2084 Int_t jVertices=0;
2085
2086 // Access to the AOD container of tracks
2087 TClonesArray &tracks = *(aod->GetTracks());
2088 Int_t jTracks=0;
2089
2090 // Add primary vertex. The primary tracks will be defined
2091 // after the loops on the composite objects (V0, cascades, kinks)
2092 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2093
2094 vtx->GetXYZ(pos); // position
2095 vtx->GetCovMatrix(covVtx); //covariance matrix
2096
2097 AliAODVertex * primary = new(vertices[jVertices++])
2098 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2099
2100 // Create vertices starting from the most complex objects
2101
2102 // Cascades
2103 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2104 AliESDcascade *cascade = esd->GetCascade(nCascade);
2105
2106 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2107 cascade->GetPosCovXi(covVtx);
2108
2109 // Add the cascade vertex
2110 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2111 covVtx,
2112 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2113 primary,
2114 AliAODVertex::kCascade);
2115
2116 primary->AddDaughter(vcascade);
2117
2118 // Add the V0 from the cascade. The ESD class have to be optimized...
2119 // Now we have to search for the corresponding Vo in the list of V0s
2120 // using the indeces of the positive and negative tracks
2121
2122 Int_t posFromV0 = cascade->GetPindex();
2123 Int_t negFromV0 = cascade->GetNindex();
2124
2125
2126 AliESDv0 * v0 = 0x0;
2127 Int_t indV0 = -1;
2128
2129 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2130
2131 v0 = esd->GetV0(iV0);
2132 Int_t posV0 = v0->GetPindex();
2133 Int_t negV0 = v0->GetNindex();
2134
2135 if (posV0==posFromV0 && negV0==negFromV0) {
2136 indV0 = iV0;
2137 break;
2138 }
2139 }
2140
2141 AliAODVertex * vV0FromCascade = 0x0;
2142
2143 if (indV0>-1 && !usedV0[indV0] ) {
2144
2145 // the V0 exists in the array of V0s and is not used
2146
2147 usedV0[indV0] = kTRUE;
2148
2149 v0->GetXYZ(pos[0], pos[1], pos[2]);
2150 v0->GetPosCov(covVtx);
2151
2152 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2153 covVtx,
2154 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2155 vcascade,
2156 AliAODVertex::kV0);
2157 } else {
2158
2159 // the V0 doesn't exist in the array of V0s or was used
2160 cerr << "Error: event " << iEvent << " cascade " << nCascade
2161 << " The V0 " << indV0
2162 << " doesn't exist in the array of V0s or was used!" << endl;
2163
2164 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2165 cascade->GetPosCov(covVtx);
2166
2167 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2168 covVtx,
2169 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2170 vcascade,
2171 AliAODVertex::kV0);
2172 vcascade->AddDaughter(vV0FromCascade);
2173 }
2174
2175 // Add the positive tracks from the V0
2176
2177 if (! usedTrack[posFromV0]) {
2178
2179 usedTrack[posFromV0] = kTRUE;
2180
2181 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2182 esdTrack->GetPxPyPz(p);
2183 esdTrack->GetXYZ(pos);
2184 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2185 esdTrack->GetESDpid(pid);
2186
2187 vV0FromCascade->AddDaughter(aodTrack =
2188 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2189 esdTrack->GetLabel(),
2190 p,
2191 kTRUE,
2192 pos,
2193 kFALSE,
2194 covTr,
2195 (Short_t)esdTrack->GetSign(),
2196 esdTrack->GetITSClusterMap(),
2197 pid,
2198 vV0FromCascade,
2199 kTRUE, // check if this is right
2200 kFALSE, // check if this is right
2201 AliAODTrack::kSecondary)
2202 );
2203 aodTrack->ConvertAliPIDtoAODPID();
2204 }
2205 else {
2206 cerr << "Error: event " << iEvent << " cascade " << nCascade
2207 << " track " << posFromV0 << " has already been used!" << endl;
2208 }
2209
2210 // Add the negative tracks from the V0
2211
2212 if (!usedTrack[negFromV0]) {
2213
2214 usedTrack[negFromV0] = kTRUE;
2215
2216 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2217 esdTrack->GetPxPyPz(p);
2218 esdTrack->GetXYZ(pos);
2219 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2220 esdTrack->GetESDpid(pid);
2221
2222 vV0FromCascade->AddDaughter(aodTrack =
2223 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2224 esdTrack->GetLabel(),
2225 p,
2226 kTRUE,
2227 pos,
2228 kFALSE,
2229 covTr,
2230 (Short_t)esdTrack->GetSign(),
2231 esdTrack->GetITSClusterMap(),
2232 pid,
2233 vV0FromCascade,
2234 kTRUE, // check if this is right
2235 kFALSE, // check if this is right
2236 AliAODTrack::kSecondary)
2237 );
2238 aodTrack->ConvertAliPIDtoAODPID();
2239 }
2240 else {
2241 cerr << "Error: event " << iEvent << " cascade " << nCascade
2242 << " track " << negFromV0 << " has already been used!" << endl;
2243 }
2244
2245 // Add the bachelor track from the cascade
2246
2247 Int_t bachelor = cascade->GetBindex();
2248
2249 if(!usedTrack[bachelor]) {
2250
2251 usedTrack[bachelor] = kTRUE;
2252
2253 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2254 esdTrack->GetPxPyPz(p);
2255 esdTrack->GetXYZ(pos);
2256 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2257 esdTrack->GetESDpid(pid);
2258
2259 vcascade->AddDaughter(aodTrack =
2260 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2261 esdTrack->GetLabel(),
2262 p,
2263 kTRUE,
2264 pos,
2265 kFALSE,
2266 covTr,
2267 (Short_t)esdTrack->GetSign(),
2268 esdTrack->GetITSClusterMap(),
2269 pid,
2270 vcascade,
2271 kTRUE, // check if this is right
2272 kFALSE, // check if this is right
2273 AliAODTrack::kSecondary)
2274 );
2275 aodTrack->ConvertAliPIDtoAODPID();
2276 }
2277 else {
2278 cerr << "Error: event " << iEvent << " cascade " << nCascade
2279 << " track " << bachelor << " has already been used!" << endl;
2280 }
2281
2282 // Add the primary track of the cascade (if any)
2283
2284 } // end of the loop on cascades
2285
2286 // V0s
2287
2288 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2289
2290 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2291
2292 AliESDv0 *v0 = esd->GetV0(nV0);
2293
2294 v0->GetXYZ(pos[0], pos[1], pos[2]);
2295 v0->GetPosCov(covVtx);
2296
2297 AliAODVertex * vV0 =
2298 new(vertices[jVertices++]) AliAODVertex(pos,
2299 covVtx,
2300 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2301 primary,
2302 AliAODVertex::kV0);
2303 primary->AddDaughter(vV0);
2304
2305 Int_t posFromV0 = v0->GetPindex();
2306 Int_t negFromV0 = v0->GetNindex();
2307
2308 // Add the positive tracks from the V0
2309
2310 if (!usedTrack[posFromV0]) {
2311
2312 usedTrack[posFromV0] = kTRUE;
2313
2314 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2315 esdTrack->GetPxPyPz(p);
2316 esdTrack->GetXYZ(pos);
2317 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2318 esdTrack->GetESDpid(pid);
2319
2320 vV0->AddDaughter(aodTrack =
2321 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2322 esdTrack->GetLabel(),
2323 p,
2324 kTRUE,
2325 pos,
2326 kFALSE,
2327 covTr,
2328 (Short_t)esdTrack->GetSign(),
2329 esdTrack->GetITSClusterMap(),
2330 pid,
2331 vV0,
2332 kTRUE, // check if this is right
2333 kFALSE, // check if this is right
2334 AliAODTrack::kSecondary)
2335 );
2336 aodTrack->ConvertAliPIDtoAODPID();
2337 }
2338 else {
2339 cerr << "Error: event " << iEvent << " V0 " << nV0
2340 << " track " << posFromV0 << " has already been used!" << endl;
2341 }
a7807689 2342
f29f1726 2343 // Add the negative tracks from the V0
2344
2345 if (!usedTrack[negFromV0]) {
2346
2347 usedTrack[negFromV0] = kTRUE;
2348
2349 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2350 esdTrack->GetPxPyPz(p);
2351 esdTrack->GetXYZ(pos);
2352 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2353 esdTrack->GetESDpid(pid);
2354
2355 vV0->AddDaughter(aodTrack =
2356 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2357 esdTrack->GetLabel(),
2358 p,
2359 kTRUE,
2360 pos,
2361 kFALSE,
2362 covTr,
2363 (Short_t)esdTrack->GetSign(),
2364 esdTrack->GetITSClusterMap(),
2365 pid,
2366 vV0,
2367 kTRUE, // check if this is right
2368 kFALSE, // check if this is right
2369 AliAODTrack::kSecondary)
2370 );
2371 aodTrack->ConvertAliPIDtoAODPID();
2372 }
2373 else {
2374 cerr << "Error: event " << iEvent << " V0 " << nV0
2375 << " track " << negFromV0 << " has already been used!" << endl;
2376 }
2377
2378 } // end of the loop on V0s
2379
2380 // Kinks: it is a big mess the access to the information in the kinks
2381 // The loop is on the tracks in order to find the mother and daugther of each kink
2382
2383
2384 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2385
2386
2387 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2388
2389 Int_t ikink = esdTrack->GetKinkIndex(0);
2390
2391 if (ikink) {
2392 // Negative kink index: mother, positive: daughter
2393
2394 // Search for the second track of the kink
2395
2396 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2397
2398 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2399
2400 Int_t jkink = esdTrack1->GetKinkIndex(0);
2401
2402 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2403
2404 // The two tracks are from the same kink
2405
2406 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2407
2408 Int_t imother = -1;
2409 Int_t idaughter = -1;
2410
2411 if (ikink<0 && jkink>0) {
2412
2413 imother = iTrack;
2414 idaughter = jTrack;
2415 }
2416 else if (ikink>0 && jkink<0) {
2417
2418 imother = jTrack;
2419 idaughter = iTrack;
2420 }
2421 else {
2422 cerr << "Error: Wrong combination of kink indexes: "
2423 << ikink << " " << jkink << endl;
2424 continue;
2425 }
2426
2427 // Add the mother track
2428
2429 AliAODTrack * mother = NULL;
2430
2431 if (!usedTrack[imother]) {
2432
2433 usedTrack[imother] = kTRUE;
2434
2435 AliESDtrack *esdTrack = esd->GetTrack(imother);
2436 esdTrack->GetPxPyPz(p);
2437 esdTrack->GetXYZ(pos);
2438 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2439 esdTrack->GetESDpid(pid);
2440
2441 mother =
2442 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2443 esdTrack->GetLabel(),
2444 p,
2445 kTRUE,
2446 pos,
2447 kFALSE,
2448 covTr,
2449 (Short_t)esdTrack->GetSign(),
2450 esdTrack->GetITSClusterMap(),
2451 pid,
2452 primary,
2453 kTRUE, // check if this is right
2454 kTRUE, // check if this is right
2455 AliAODTrack::kPrimary);
2456 primary->AddDaughter(mother);
2457 mother->ConvertAliPIDtoAODPID();
2458 }
2459 else {
2460 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2461 << " track " << imother << " has already been used!" << endl;
2462 }
2463
2464 // Add the kink vertex
2465 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2466
2467 AliAODVertex * vkink =
2468 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2469 NULL,
2470 0.,
2471 mother,
2472 AliAODVertex::kKink);
2473 // Add the daughter track
2474
2475 AliAODTrack * daughter = NULL;
2476
2477 if (!usedTrack[idaughter]) {
2478
2479 usedTrack[idaughter] = kTRUE;
2480
2481 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2482 esdTrack->GetPxPyPz(p);
2483 esdTrack->GetXYZ(pos);
2484 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2485 esdTrack->GetESDpid(pid);
2486
2487 daughter =
2488 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2489 esdTrack->GetLabel(),
2490 p,
2491 kTRUE,
2492 pos,
2493 kFALSE,
2494 covTr,
2495 (Short_t)esdTrack->GetSign(),
2496 esdTrack->GetITSClusterMap(),
2497 pid,
2498 vkink,
2499 kTRUE, // check if this is right
2500 kTRUE, // check if this is right
2501 AliAODTrack::kPrimary);
2502 vkink->AddDaughter(daughter);
2503 daughter->ConvertAliPIDtoAODPID();
2504 }
2505 else {
2506 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2507 << " track " << idaughter << " has already been used!" << endl;
2508 }
2509
2510
2511 }
2512 }
2513
2514 }
2515
2516 }
2517
2518
2519 // Tracks (primary and orphan)
2520
2521 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2522
2523
2524 if (usedTrack[nTrack]) continue;
2525
2526 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2527 esdTrack->GetPxPyPz(p);
2528 esdTrack->GetXYZ(pos);
2529 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2530 esdTrack->GetESDpid(pid);
2531
2532 Float_t impactXY, impactZ;
2533
2534 esdTrack->GetImpactParameters(impactXY,impactZ);
2535
2536 if (impactXY<3) {
2537 // track inside the beam pipe
2538
2539 primary->AddDaughter(aodTrack =
2540 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2541 esdTrack->GetLabel(),
2542 p,
2543 kTRUE,
2544 pos,
2545 kFALSE,
2546 covTr,
2547 (Short_t)esdTrack->GetSign(),
2548 esdTrack->GetITSClusterMap(),
2549 pid,
2550 primary,
2551 kTRUE, // check if this is right
2552 kTRUE, // check if this is right
2553 AliAODTrack::kPrimary)
2554 );
2555 aodTrack->ConvertAliPIDtoAODPID();
2556 }
2557 else {
2558 // outside the beam pipe: orphan track
2559 aodTrack =
2560 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2561 esdTrack->GetLabel(),
2562 p,
2563 kTRUE,
2564 pos,
2565 kFALSE,
2566 covTr,
2567 (Short_t)esdTrack->GetSign(),
2568 esdTrack->GetITSClusterMap(),
2569 pid,
2570 NULL,
2571 kFALSE, // check if this is right
2572 kFALSE, // check if this is right
2573 AliAODTrack::kOrphan);
2574 aodTrack->ConvertAliPIDtoAODPID();
2575 }
2576 } // end of loop on tracks
2577
2578 // muon tracks
2579 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2580 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2581
2582 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2583 p[0] = esdMuTrack->Px();
2584 p[1] = esdMuTrack->Py();
2585 p[2] = esdMuTrack->Pz();
2586 pos[0] = primary->GetX();
2587 pos[1] = primary->GetY();
2588 pos[2] = primary->GetZ();
2589
2590 // has to be changed once the muon pid is provided by the ESD
2591 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2592
2593 primary->AddDaughter(
2594 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2595 0, // no label provided
2596 p,
2597 kTRUE,
2598 pos,
2599 kFALSE,
2600 NULL, // no covariance matrix provided
2601 (Short_t)-99, // no charge provided
2602 0, // no ITSClusterMap
2603 pid,
2604 primary,
2605 kTRUE, // check if this is right
2606 kTRUE, // not used for vertex fit
2607 AliAODTrack::kPrimary)
2608 );
2609 }
2610
2611 // Access to the AOD container of clusters
2612 TClonesArray &clusters = *(aod->GetClusters());
2613 Int_t jClusters=0;
2614
2615 // Calo Clusters
2616 Int_t nClusters = esd->GetNumberOfCaloClusters();
2617
2618 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2619
2620 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2621
2622 Int_t id = cluster->GetID();
2623 Int_t label = -1;
2624 Float_t energy = cluster->GetClusterEnergy();
2625 cluster->GetGlobalPosition(posF);
2626 AliAODVertex *prodVertex = primary;
2627 AliAODTrack *primTrack = NULL;
2628 Char_t ttype=AliAODCluster::kUndef;
2629
2630 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2631 else if (cluster->IsEMCAL()) {
2632
2633 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2634 ttype = AliAODCluster::kEMCALPseudoCluster;
2635 else
2636 ttype = AliAODCluster::kEMCALClusterv1;
2637
2638 }
2639
2640 new(clusters[jClusters++]) AliAODCluster(id,
2641 label,
2642 energy,
2643 pos,
2644 NULL, // no covariance matrix provided
2645 NULL, // no pid for clusters provided
2646 prodVertex,
2647 primTrack,
2648 ttype);
2649
2650 } // end of loop on calo clusters
2651
2652 delete [] usedTrack;
2653 delete [] usedV0;
2654 delete [] usedKink;
2655
2656 // fill the tree for this event
2657 aodTree->Fill();
2658 } // end of event loop
2659
2660 aodTree->GetUserInfo()->Add(aod);
2661
2662 // close ESD file
2663 esdFile->Close();
2664
2665 // write the tree to the specified file
2666 aodFile = aodTree->GetCurrentFile();
2667 aodFile->cd();
2668 aodTree->Write();
2669
a7807689 2670 return;
2671}
2672
2673
98937d93 2674void AliReconstruction::WriteAlignmentData(AliESD* esd)
2675{
2676 // Write space-points which are then used in the alignment procedures
2677 // For the moment only ITS, TRD and TPC
2678
2679 // Load TOF clusters
d528ee75 2680 if (fTracker[3]){
2681 fLoader[3]->LoadRecPoints("read");
2682 TTree* tree = fLoader[3]->TreeR();
2683 if (!tree) {
2684 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2685 return;
2686 }
2687 fTracker[3]->LoadClusters(tree);
98937d93 2688 }
98937d93 2689 Int_t ntracks = esd->GetNumberOfTracks();
2690 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2691 {
2692 AliESDtrack *track = esd->GetTrack(itrack);
2693 Int_t nsp = 0;
ef7253ac 2694 Int_t idx[200];
98937d93 2695 for (Int_t iDet = 3; iDet >= 0; iDet--)
2696 nsp += track->GetNcls(iDet);
2697 if (nsp) {
2698 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2699 track->SetTrackPointArray(sp);
2700 Int_t isptrack = 0;
2701 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2702 AliTracker *tracker = fTracker[iDet];
2703 if (!tracker) continue;
2704 Int_t nspdet = track->GetNcls(iDet);
98937d93 2705 if (nspdet <= 0) continue;
2706 track->GetClusters(iDet,idx);
2707 AliTrackPoint p;
2708 Int_t isp = 0;
2709 Int_t isp2 = 0;
2710 while (isp < nspdet) {
2711 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 2712 const Int_t kNTPCmax = 159;
2713 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2714 if (!isvalid) continue;
2715 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2716 }
98937d93 2717 }
2718 }
2719 }
d528ee75 2720 if (fTracker[3]){
2721 fTracker[3]->UnloadClusters();
2722 fLoader[3]->UnloadRecPoints();
2723 }
98937d93 2724}
2e3550da 2725
2726//_____________________________________________________________________________
2727void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2728{
2729 // The method reads the raw-data error log
2730 // accumulated within the rawReader.
2731 // It extracts the raw-data errors related to
2732 // the current event and stores them into
2733 // a TClonesArray inside the esd object.
2734
2735 if (!fRawReader) return;
2736
2737 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2738
2739 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2740 if (!log) continue;
2741 if (iEvent != log->GetEventNumber()) continue;
2742
2743 esd->AddRawDataErrorLog(log);
2744 }
2745
2746}