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