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