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