Separate methods Init(Option_t *) (called by AnalysisManager::SlaveBegin)
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
... / ...
CommitLineData
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. //
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// //
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// //
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// //
59// The name of the galice file can be changed from the default //
60// "galice.root" by passing it as argument to the AliReconstruction //
61// constructor or by //
62// //
63// rec.SetGAliceFile("..."); //
64// //
65// The local reconstruction can be switched on or off for individual //
66// detectors by //
67// //
68// rec.SetRunLocalReconstruction("..."); //
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// //
74// The reconstruction of the primary vertex position can be switched off by //
75// //
76// rec.SetRunVertexFinder(kFALSE); //
77// //
78// The tracking and the creation of ESD tracks can be switched on for //
79// selected detectors by //
80// //
81// rec.SetRunTracking("..."); //
82// //
83// Uniform/nonuniform field tracking switches (default: uniform field) //
84// //
85// rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
86// //
87// The filling of additional ESD information can be steered by //
88// //
89// rec.SetFillESD("..."); //
90// //
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. //
100// //
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. //
103// //
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// //
115///////////////////////////////////////////////////////////////////////////////
116
117#include <TArrayF.h>
118#include <TFile.h>
119#include <TList.h>
120#include <TSystem.h>
121#include <TROOT.h>
122#include <TPluginManager.h>
123#include <TGeoManager.h>
124#include <TLorentzVector.h>
125#include <TArrayS.h>
126#include <TArrayD.h>
127
128#include "AliReconstruction.h"
129#include "AliCodeTimer.h"
130#include "AliReconstructor.h"
131#include "AliLog.h"
132#include "AliRunLoader.h"
133#include "AliRun.h"
134#include "AliRawReaderFile.h"
135#include "AliRawReaderDate.h"
136#include "AliRawReaderRoot.h"
137#include "AliRawEventHeaderBase.h"
138#include "AliESDEvent.h"
139#include "AliESDMuonTrack.h"
140#include "AliESDfriend.h"
141#include "AliESDVertex.h"
142#include "AliESDcascade.h"
143#include "AliESDkink.h"
144#include "AliESDtrack.h"
145#include "AliESDCaloCluster.h"
146#include "AliMultiplicity.h"
147#include "AliTracker.h"
148#include "AliVertexer.h"
149#include "AliVertexerTracks.h"
150#include "AliV0vertexer.h"
151#include "AliCascadeVertexer.h"
152#include "AliHeader.h"
153#include "AliGenEventHeader.h"
154#include "AliPID.h"
155#include "AliESDpid.h"
156#include "AliESDtrack.h"
157#include "AliESDPmdTrack.h"
158
159#include "AliESDTagCreator.h"
160#include "AliAODTagCreator.h"
161
162#include "AliGeomManager.h"
163#include "AliTrackPointArray.h"
164#include "AliCDBManager.h"
165#include "AliCDBStorage.h"
166#include "AliCDBEntry.h"
167#include "AliAlignObj.h"
168
169#include "AliCentralTrigger.h"
170#include "AliCTPRawStream.h"
171
172#include "AliAODEvent.h"
173#include "AliAODHeader.h"
174#include "AliAODTrack.h"
175#include "AliAODVertex.h"
176#include "AliAODv0.h"
177#include "AliAODJet.h"
178#include "AliAODCaloCells.h"
179#include "AliAODCaloCluster.h"
180#include "AliAODPmdCluster.h"
181#include "AliAODFmdCluster.h"
182#include "AliAODTracklets.h"
183
184#include "AliQADataMaker.h"
185#include "AliQA.h"
186#include "AliQADataMakerSteer.h"
187
188#include "AliSysInfo.h" // memory snapshots
189
190
191ClassImp(AliReconstruction)
192
193
194//_____________________________________________________________________________
195const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
196
197//_____________________________________________________________________________
198AliReconstruction::AliReconstruction(const char* gAliceFilename,
199 const char* name, const char* title) :
200 TNamed(name, title),
201
202 fUniformField(kTRUE),
203 fRunVertexFinder(kTRUE),
204 fRunHLTTracking(kFALSE),
205 fRunMuonTracking(kFALSE),
206 fRunV0Finder(kTRUE),
207 fRunCascadeFinder(kTRUE),
208 fStopOnError(kFALSE),
209 fWriteAlignmentData(kFALSE),
210 fWriteESDfriend(kFALSE),
211 fWriteAOD(kFALSE),
212 fFillTriggerESD(kTRUE),
213
214 fCleanESD(kTRUE),
215 fV0DCAmax(3.),
216 fV0CsPmin(0.),
217 fDmax(50.),
218 fZmax(50.),
219
220 fRunLocalReconstruction("ALL"),
221 fRunTracking("ALL"),
222 fFillESD("ALL"),
223 fUseTrackingErrorsForAlignment(""),
224 fGAliceFileName(gAliceFilename),
225 fInput(""),
226 fEquipIdMap(""),
227 fFirstEvent(0),
228 fLastEvent(-1),
229 fNumberOfEventsPerFile(1),
230 fCheckPointLevel(0),
231 fOptions(),
232 fLoadAlignFromCDB(kTRUE),
233 fLoadAlignData("ALL"),
234 fESDPar(""),
235
236 fRunLoader(NULL),
237 fRawReader(NULL),
238
239 fVertexer(NULL),
240 fDiamondProfile(NULL),
241
242 fGRPList(NULL),
243
244 fAlignObjArray(NULL),
245 fCDBUri(),
246 fSpecCDBUri(),
247 fInitCDBCalled(kFALSE),
248 fSetRunNumberFromDataCalled(kFALSE),
249 fRunQA(kTRUE),
250 fQAInLoop(kFALSE)
251
252{
253// create reconstruction object with default parameters
254
255 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
256 fReconstructor[iDet] = NULL;
257 fLoader[iDet] = NULL;
258 fTracker[iDet] = NULL;
259 fQADataMaker[iDet] = NULL;
260 fQACycles[iDet] = 999999;
261 }
262 AliPID pid;
263}
264
265//_____________________________________________________________________________
266AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
267 TNamed(rec),
268
269 fUniformField(rec.fUniformField),
270 fRunVertexFinder(rec.fRunVertexFinder),
271 fRunHLTTracking(rec.fRunHLTTracking),
272 fRunMuonTracking(rec.fRunMuonTracking),
273 fRunV0Finder(rec.fRunV0Finder),
274 fRunCascadeFinder(rec.fRunCascadeFinder),
275 fStopOnError(rec.fStopOnError),
276 fWriteAlignmentData(rec.fWriteAlignmentData),
277 fWriteESDfriend(rec.fWriteESDfriend),
278 fWriteAOD(rec.fWriteAOD),
279 fFillTriggerESD(rec.fFillTriggerESD),
280
281 fCleanESD(rec.fCleanESD),
282 fV0DCAmax(rec.fV0DCAmax),
283 fV0CsPmin(rec.fV0CsPmin),
284 fDmax(rec.fDmax),
285 fZmax(rec.fZmax),
286
287 fRunLocalReconstruction(rec.fRunLocalReconstruction),
288 fRunTracking(rec.fRunTracking),
289 fFillESD(rec.fFillESD),
290 fUseTrackingErrorsForAlignment(rec.fUseTrackingErrorsForAlignment),
291 fGAliceFileName(rec.fGAliceFileName),
292 fInput(rec.fInput),
293 fEquipIdMap(rec.fEquipIdMap),
294 fFirstEvent(rec.fFirstEvent),
295 fLastEvent(rec.fLastEvent),
296 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
297 fCheckPointLevel(0),
298 fOptions(),
299 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
300 fLoadAlignData(rec.fLoadAlignData),
301 fESDPar(rec.fESDPar),
302
303 fRunLoader(NULL),
304 fRawReader(NULL),
305
306 fVertexer(NULL),
307 fDiamondProfile(NULL),
308
309 fGRPList(NULL),
310
311 fAlignObjArray(rec.fAlignObjArray),
312 fCDBUri(rec.fCDBUri),
313 fSpecCDBUri(),
314 fInitCDBCalled(rec.fInitCDBCalled),
315 fSetRunNumberFromDataCalled(rec.fSetRunNumberFromDataCalled),
316 fRunQA(kTRUE),
317 fQAInLoop(kFALSE)
318
319{
320// copy constructor
321
322 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
323 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
324 }
325 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
326 fReconstructor[iDet] = NULL;
327 fLoader[iDet] = NULL;
328 fTracker[iDet] = NULL;
329 fQADataMaker[iDet] = NULL;
330 fQACycles[iDet] = rec.fQACycles[iDet];
331 }
332 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
333 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
334 }
335}
336
337//_____________________________________________________________________________
338AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
339{
340// assignment operator
341
342 this->~AliReconstruction();
343 new(this) AliReconstruction(rec);
344 return *this;
345}
346
347//_____________________________________________________________________________
348AliReconstruction::~AliReconstruction()
349{
350// clean up
351
352 CleanUp();
353 fOptions.Delete();
354 fSpecCDBUri.Delete();
355
356 AliCodeTimer::Instance()->Print();
357}
358
359//_____________________________________________________________________________
360void AliReconstruction::InitCDB()
361{
362// activate a default CDB storage
363// First check if we have any CDB storage set, because it is used
364// to retrieve the calibration and alignment constants
365
366 if (fInitCDBCalled) return;
367 fInitCDBCalled = kTRUE;
368
369 AliCDBManager* man = AliCDBManager::Instance();
370 if (man->IsDefaultStorageSet())
371 {
372 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
373 AliWarning("Default CDB storage has been already set !");
374 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
375 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
376 fCDBUri = man->GetDefaultStorage()->GetURI();
377 }
378 else {
379 if (fCDBUri.Length() > 0)
380 {
381 AliDebug(2,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
382 AliDebug(2, Form("Default CDB storage is set to: %s", fCDBUri.Data()));
383 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
384 } else {
385 fCDBUri="local://$ALICE_ROOT";
386 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
387 AliWarning("Default CDB storage not yet set !!!!");
388 AliWarning(Form("Setting it now to: %s", fCDBUri.Data()));
389 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
390
391 }
392 man->SetDefaultStorage(fCDBUri);
393 }
394
395 // Now activate the detector specific CDB storage locations
396 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
397 TObject* obj = fSpecCDBUri[i];
398 if (!obj) continue;
399 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
400 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
401 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
402 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
403 }
404
405}
406
407//_____________________________________________________________________________
408void AliReconstruction::SetDefaultStorage(const char* uri) {
409// Store the desired default CDB storage location
410// Activate it later within the Run() method
411
412 fCDBUri = uri;
413
414}
415
416//_____________________________________________________________________________
417void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
418// Store a detector-specific CDB storage location
419// Activate it later within the Run() method
420
421 AliCDBPath aPath(calibType);
422 if(!aPath.IsValid()){
423 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
424 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
425 if(!strcmp(calibType, fgkDetectorName[iDet])) {
426 aPath.SetPath(Form("%s/*", calibType));
427 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
428 break;
429 }
430 }
431 if(!aPath.IsValid()){
432 AliError(Form("Not a valid path or detector: %s", calibType));
433 return;
434 }
435 }
436
437// // check that calibType refers to a "valid" detector name
438// Bool_t isDetector = kFALSE;
439// for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
440// TString detName = fgkDetectorName[iDet];
441// if(aPath.GetLevel0() == detName) {
442// isDetector = kTRUE;
443// break;
444// }
445// }
446//
447// if(!isDetector) {
448// AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
449// return;
450// }
451
452 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
453 if (obj) fSpecCDBUri.Remove(obj);
454 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
455
456}
457
458//_____________________________________________________________________________
459Bool_t AliReconstruction::SetRunNumberFromData()
460{
461 // The method is called in Run() in order
462 // to set a correct run number.
463 // In case of raw data reconstruction the
464 // run number is taken from the raw data header
465
466 if (fSetRunNumberFromDataCalled) return kTRUE;
467 fSetRunNumberFromDataCalled = kTRUE;
468
469 AliCDBManager* man = AliCDBManager::Instance();
470
471 if(man->GetRun() > 0) {
472 AliWarning("Run number is taken from event header! Ignoring settings in AliCDBManager!");
473 }
474
475 if (!fRunLoader) {
476 AliError("No run loader is found !");
477 return kFALSE;
478 }
479 // read run number from gAlice
480 if(fRunLoader->GetAliRun())
481 AliCDBManager::Instance()->SetRun(fRunLoader->GetHeader()->GetRun());
482 else {
483 if(fRawReader) {
484 if(fRawReader->NextEvent()) {
485 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
486 fRawReader->RewindEvents();
487 }
488 else {
489 AliError("No raw-data events found !");
490 return kFALSE;
491 }
492 }
493 else {
494 AliError("Neither gAlice nor RawReader objects are found !");
495 return kFALSE;
496 }
497 }
498
499 man->Print();
500
501 return kTRUE;
502}
503
504//_____________________________________________________________________________
505void AliReconstruction::SetCDBLock() {
506 // Set CDB lock: from now on it is forbidden to reset the run number
507 // or the default storage or to activate any further storage!
508
509 AliCDBManager::Instance()->SetLock(1);
510}
511
512//_____________________________________________________________________________
513Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
514{
515 // Read the alignment objects from CDB.
516 // Each detector is supposed to have the
517 // alignment objects in DET/Align/Data CDB path.
518 // All the detector objects are then collected,
519 // sorted by geometry level (starting from ALIC) and
520 // then applied to the TGeo geometry.
521 // Finally an overlaps check is performed.
522
523 // Load alignment data from CDB and fill fAlignObjArray
524 if(fLoadAlignFromCDB){
525
526 TString detStr = detectors;
527 TString loadAlObjsListOfDets = "";
528
529 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
530 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
531 loadAlObjsListOfDets += fgkDetectorName[iDet];
532 loadAlObjsListOfDets += " ";
533 } // end loop over detectors
534 loadAlObjsListOfDets.Prepend("GRP "); //add alignment objects for non-sensitive modules
535 AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
536 }else{
537 // Check if the array with alignment objects was
538 // provided by the user. If yes, apply the objects
539 // to the present TGeo geometry
540 if (fAlignObjArray) {
541 if (gGeoManager && gGeoManager->IsClosed()) {
542 if (AliGeomManager::ApplyAlignObjsToGeom(*fAlignObjArray) == kFALSE) {
543 AliError("The misalignment of one or more volumes failed!"
544 "Compare the list of simulated detectors and the list of detector alignment data!");
545 return kFALSE;
546 }
547 }
548 else {
549 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
550 return kFALSE;
551 }
552 }
553 }
554
555 delete fAlignObjArray; fAlignObjArray=0;
556
557 return kTRUE;
558}
559
560//_____________________________________________________________________________
561void AliReconstruction::SetGAliceFile(const char* fileName)
562{
563// set the name of the galice file
564
565 fGAliceFileName = fileName;
566}
567
568//_____________________________________________________________________________
569void AliReconstruction::SetOption(const char* detector, const char* option)
570{
571// set options for the reconstruction of a detector
572
573 TObject* obj = fOptions.FindObject(detector);
574 if (obj) fOptions.Remove(obj);
575 fOptions.Add(new TNamed(detector, option));
576}
577
578
579//_____________________________________________________________________________
580Bool_t AliReconstruction::Run(const char* input)
581{
582// run the reconstruction
583
584 AliCodeTimerAuto("")
585
586 // set the input
587 if (!input) input = fInput.Data();
588 TString fileName(input);
589 if (fileName.EndsWith("/")) {
590 fRawReader = new AliRawReaderFile(fileName);
591 } else if (fileName.EndsWith(".root")) {
592 fRawReader = new AliRawReaderRoot(fileName);
593 } else if (!fileName.IsNull()) {
594 fRawReader = new AliRawReaderDate(fileName);
595 fRawReader->SelectEvents(7);
596 }
597 if (!fEquipIdMap.IsNull() && fRawReader)
598 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
599
600 AliSysInfo::AddStamp("Start");
601 // get the run loader
602 if (!InitRunLoader()) return kFALSE;
603 AliSysInfo::AddStamp("LoadLoader");
604
605 // Initialize the CDB storage
606 InitCDB();
607
608 AliSysInfo::AddStamp("LoadCDB");
609
610 // Set run number in CDBManager (if it is not already set by the user)
611 if (!SetRunNumberFromData()) if (fStopOnError) return kFALSE;
612
613 // Set CDB lock: from now on it is forbidden to reset the run number
614 // or the default storage or to activate any further storage!
615 SetCDBLock();
616
617 // Import ideal TGeo geometry and apply misalignment
618 if (!gGeoManager) {
619 TString geom(gSystem->DirName(fGAliceFileName));
620 geom += "/geometry.root";
621 AliGeomManager::LoadGeometry(geom.Data());
622 if (!gGeoManager) if (fStopOnError) return kFALSE;
623 }
624
625 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
626 AliSysInfo::AddStamp("LoadGeom");
627
628 //QA
629 AliQADataMakerSteer qas ;
630 if ( fRunQA && fRawReader)
631 qas.Run("ALL", fRawReader) ;
632
633 // checking the QA of previous steps
634 CheckQA() ;
635
636 // local reconstruction
637 if (!fRunLocalReconstruction.IsNull()) {
638 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
639 if (fStopOnError) {CleanUp(); return kFALSE;}
640 }
641 }
642// if (!fRunVertexFinder && fRunTracking.IsNull() &&
643// fFillESD.IsNull()) return kTRUE;
644
645 // get vertexer
646 if (fRunVertexFinder && !CreateVertexer()) {
647 if (fStopOnError) {
648 CleanUp();
649 return kFALSE;
650 }
651 }
652 AliSysInfo::AddStamp("Vertexer");
653
654 // get trackers
655 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
656 if (fStopOnError) {
657 CleanUp();
658 return kFALSE;
659 }
660 }
661 AliSysInfo::AddStamp("LoadTrackers");
662
663 // get the possibly already existing ESD file and tree
664 AliESDEvent* esd = new AliESDEvent(); AliESDEvent* hltesd = new AliESDEvent();
665 TFile* fileOld = NULL;
666 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
667 if (!gSystem->AccessPathName("AliESDs.root")){
668 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
669 fileOld = TFile::Open("AliESDs.old.root");
670 if (fileOld && fileOld->IsOpen()) {
671 treeOld = (TTree*) fileOld->Get("esdTree");
672 if (treeOld)esd->ReadFromTree(treeOld);
673 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
674 if (hlttreeOld) hltesd->ReadFromTree(hlttreeOld);
675 }
676 }
677
678 // create the ESD output file and tree
679 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
680 file->SetCompressionLevel(2);
681 if (!file->IsOpen()) {
682 AliError("opening AliESDs.root failed");
683 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
684 }
685
686 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
687 esd = new AliESDEvent();
688 esd->CreateStdContent();
689 esd->WriteToTree(tree);
690
691 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
692 hltesd = new AliESDEvent();
693 hltesd->CreateStdContent();
694 hltesd->WriteToTree(hlttree);
695
696 /* CKB Why?
697 delete esd; delete hltesd;
698 esd = NULL; hltesd = NULL;
699 */
700 // create the branch with ESD additions
701
702
703
704 AliESDfriend *esdf = 0;
705 if (fWriteESDfriend) {
706 esdf = new AliESDfriend();
707 TBranch *br=tree->Branch("ESDfriend.","AliESDfriend", &esdf);
708 br->SetFile("AliESDfriends.root");
709 esd->AddObject(esdf);
710 }
711
712
713 // Get the GRP CDB entry
714 AliCDBEntry* entryGRP = AliCDBManager::Instance()->Get("GRP/GRP/Data");
715
716 if(entryGRP) {
717 fGRPList = dynamic_cast<TList*> (entryGRP->GetObject());
718 } else {
719 AliError("No GRP entry found in OCDB!");
720 }
721
722 // Get the diamond profile from OCDB
723 AliCDBEntry* entry = AliCDBManager::Instance()
724 ->Get("GRP/Calib/MeanVertex");
725
726 if(entry) {
727 fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
728 } else {
729 AliError("No diamond profile found in OCDB!");
730 }
731
732 AliVertexerTracks tVertexer(AliTracker::GetBz());
733 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
734
735 // loop over events
736
737 if (fRawReader) fRawReader->RewindEvents();
738 TString detStr(fFillESD) ;
739
740 ProcInfo_t ProcInfo;
741 gSystem->GetProcInfo(&ProcInfo);
742 AliInfo(Form("Current memory usage %d %d", ProcInfo.fMemResident, ProcInfo.fMemVirtual));
743
744 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
745 if (fRawReader) fRawReader->NextEvent();
746 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
747 // copy old ESD to the new one
748 if (treeOld) {
749 esd->ReadFromTree(treeOld);
750 treeOld->GetEntry(iEvent);
751 }
752 tree->Fill();
753 if (hlttreeOld) {
754 esd->ReadFromTree(hlttreeOld);
755 hlttreeOld->GetEntry(iEvent);
756 }
757 hlttree->Fill();
758 continue;
759 }
760
761 AliInfo(Form("processing event %d", iEvent));
762 fRunLoader->GetEvent(iEvent);
763
764 char aFileName[256];
765 sprintf(aFileName, "ESD_%d.%d_final.root",
766 fRunLoader->GetHeader()->GetRun(),
767 fRunLoader->GetHeader()->GetEventNrInRun());
768 if (!gSystem->AccessPathName(aFileName)) continue;
769
770 // local reconstruction
771 if (!fRunLocalReconstruction.IsNull()) {
772 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
773 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
774 }
775 }
776
777
778 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
779 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
780 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
781 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
782
783 // Set magnetic field from the tracker
784 esd->SetMagneticField(AliTracker::GetBz());
785 hltesd->SetMagneticField(AliTracker::GetBz());
786
787
788
789 // Fill raw-data error log into the ESD
790 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
791
792 // vertex finder
793 if (fRunVertexFinder) {
794 if (!ReadESD(esd, "vertex")) {
795 if (!RunVertexFinder(esd)) {
796 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
797 }
798 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
799 }
800 }
801
802 // HLT tracking
803 if (!fRunTracking.IsNull()) {
804 if (fRunHLTTracking) {
805 hltesd->SetVertex(esd->GetVertex());
806 if (!RunHLTTracking(hltesd)) {
807 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
808 }
809 }
810 }
811
812 // Muon tracking
813 if (!fRunTracking.IsNull()) {
814 if (fRunMuonTracking) {
815 if (!RunMuonTracking(esd)) {
816 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
817 }
818 }
819 }
820
821 // barrel tracking
822 if (!fRunTracking.IsNull()) {
823 if (!ReadESD(esd, "tracking")) {
824 if (!RunTracking(esd)) {
825 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
826 }
827 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
828 }
829 }
830
831 // fill ESD
832 if (!fFillESD.IsNull()) {
833 if (!FillESD(esd, fFillESD)) {
834 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
835 }
836 }
837
838// if (!fFillESD.IsNull())
839
840 // do QA in the event loop if requested
841 if (fQAInLoop)
842 RunQA(fFillESD.Data(), esd);
843
844 // fill Event header information from the RawEventHeader
845 if (fRawReader){FillRawEventHeaderESD(esd);}
846
847 // combined PID
848 AliESDpid::MakePID(esd);
849 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
850
851 if (fFillTriggerESD) {
852 if (!ReadESD(esd, "trigger")) {
853 if (!FillTriggerESD(esd)) {
854 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
855 }
856 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
857 }
858 }
859
860 file->cd();
861
862 //Try to improve the reconstructed primary vertex position using the tracks
863 AliESDVertex *pvtx=0;
864 Bool_t dovertex=kTRUE;
865 TObject* obj = fOptions.FindObject("ITS");
866 if (obj) {
867 TString optITS = obj->GetTitle();
868 if (optITS.Contains("cosmics") || optITS.Contains("COSMICS"))
869 dovertex=kFALSE;
870 }
871 if(dovertex) pvtx=tVertexer.FindPrimaryVertex(esd);
872 if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
873
874 if (pvtx)
875 if (pvtx->GetStatus()) {
876 // Store the improved primary vertex
877 esd->SetPrimaryVertex(pvtx);
878 // Propagate the tracks to the DCA to the improved primary vertex
879 Double_t somethingbig = 777.;
880 Double_t bz = esd->GetMagneticField();
881 Int_t nt=esd->GetNumberOfTracks();
882 while (nt--) {
883 AliESDtrack *t = esd->GetTrack(nt);
884 t->RelateToVertex(pvtx, bz, somethingbig);
885 }
886 }
887
888 if (fRunV0Finder) {
889 // V0 finding
890 AliV0vertexer vtxer;
891 vtxer.Tracks2V0vertices(esd);
892
893 if (fRunCascadeFinder) {
894 // Cascade finding
895 AliCascadeVertexer cvtxer;
896 cvtxer.V0sTracks2CascadeVertices(esd);
897 }
898 }
899
900 // write ESD
901 if (fCleanESD) CleanESD(esd);
902 if (fWriteESDfriend) {
903 esdf->~AliESDfriend();
904 new (esdf) AliESDfriend(); // Reset...
905 esd->GetESDfriend(esdf);
906 }
907 tree->Fill();
908
909 // write HLT ESD
910 hlttree->Fill();
911
912 if (fCheckPointLevel > 0) WriteESD(esd, "final");
913 esd->Reset();
914 hltesd->Reset();
915 if (fWriteESDfriend) {
916 esdf->~AliESDfriend();
917 new (esdf) AliESDfriend(); // Reset...
918 }
919 // esdf->Reset();
920 // delete esdf; esdf = 0;
921 // ESD QA
922
923 gSystem->GetProcInfo(&ProcInfo);
924 AliInfo(Form("Event %d -> Current memory usage %d %d",iEvent, ProcInfo.fMemResident, ProcInfo.fMemVirtual));
925 }
926
927 detStr = fFillESD ;
928
929
930 // write quality assurance ESDs data (one entry for all events)
931 if (fQAInLoop) {
932 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
933 if (!IsSelected(fgkDetectorName[iDet], detStr))
934 continue;
935 AliQADataMaker * qadm = GetQADataMaker(iDet);
936 if (!qadm) continue;
937 qadm->EndOfCycle(AliQA::kRECPOINTS);
938 qadm->EndOfCycle(AliQA::kESDS);
939 qadm->Finish();
940 }
941 }
942
943 tree->GetUserInfo()->Add(esd);
944 hlttree->GetUserInfo()->Add(hltesd);
945
946 const TMap *cdbMap = AliCDBManager::Instance()->GetStorageMap();
947 const TList *cdbList = AliCDBManager::Instance()->GetRetrievedIds();
948
949 TMap *cdbMapCopy = new TMap(cdbMap->GetEntries());
950 cdbMapCopy->SetOwner(1);
951 cdbMapCopy->SetName("cdbMap");
952 TIter iter(cdbMap->GetTable());
953
954 TPair* pair = 0;
955 while((pair = dynamic_cast<TPair*> (iter.Next()))){
956 TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
957 TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
958 cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
959 }
960
961 TList *cdbListCopy = new TList();
962 cdbListCopy->SetOwner(1);
963 cdbListCopy->SetName("cdbList");
964
965 TIter iter2(cdbList);
966
967 AliCDBId* id=0;
968 while((id = dynamic_cast<AliCDBId*> (iter2.Next()))){
969 cdbListCopy->Add(new TObjString(id->ToString().Data()));
970 }
971
972 tree->GetUserInfo()->Add(cdbMapCopy);
973 tree->GetUserInfo()->Add(cdbListCopy);
974
975
976 if(fESDPar.Contains("ESD.par")){
977 AliInfo("Attaching ESD.par to Tree");
978 TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
979 tree->GetUserInfo()->Add(fn);
980 }
981
982
983 file->cd();
984 if (fWriteESDfriend)
985 tree->SetBranchStatus("ESDfriend*",0);
986 // we want to have only one tree version number
987 tree->Write(tree->GetName(),TObject::kOverwrite);
988 hlttree->Write();
989
990 if (fWriteAOD) {
991 TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
992 ESDFile2AODFile(file, aodFile);
993 aodFile->Close();
994 }
995
996 gROOT->cd();
997 CleanUp(file, fileOld);
998
999 // Create tags for the events in the ESD tree (the ESD tree is always present)
1000 // In case of empty events the tags will contain dummy values
1001 AliESDTagCreator *esdtagCreator = new AliESDTagCreator();
1002 esdtagCreator->CreateESDTags(fFirstEvent,fLastEvent,fGRPList);
1003 if (fWriteAOD) {
1004 AliAODTagCreator *aodtagCreator = new AliAODTagCreator();
1005 aodtagCreator->CreateAODTags(fFirstEvent,fLastEvent,fGRPList);
1006 }
1007
1008 //QA fone outside the event loop
1009 if (!fQAInLoop) {
1010 if ( fRunQA ) {
1011 qas.Run(fRunLocalReconstruction.Data(), AliQA::kRECPOINTS) ;
1012 // qas.Reset() ;
1013 qas.Run(fRunTracking.Data(), AliQA::kESDS) ;
1014 }
1015 }
1016
1017 // Cleanup of CDB manager: cache and active storages!
1018 AliCDBManager::Instance()->ClearCache();
1019
1020
1021 return kTRUE;
1022}
1023
1024
1025//_____________________________________________________________________________
1026Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
1027{
1028// run the local reconstruction
1029 static Int_t eventNr=0;
1030 AliCodeTimerAuto("")
1031
1032 // AliCDBManager* man = AliCDBManager::Instance();
1033// Bool_t origCache = man->GetCacheFlag();
1034
1035// TString detStr = detectors;
1036// for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1037// if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1038// AliReconstructor* reconstructor = GetReconstructor(iDet);
1039// if (!reconstructor) continue;
1040// if (reconstructor->HasLocalReconstruction()) continue;
1041
1042// AliCodeTimerStart(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1043// AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1044
1045// AliCodeTimerStart(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
1046// AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
1047
1048// man->SetCacheFlag(kTRUE);
1049// TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
1050// man->GetAll(calibPath); // entries are cached!
1051
1052// AliCodeTimerStop(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
1053
1054// if (fRawReader) {
1055// fRawReader->RewindEvents();
1056// reconstructor->Reconstruct(fRunLoader, fRawReader);
1057// } else {
1058// reconstructor->Reconstruct(fRunLoader);
1059// }
1060
1061// AliCodeTimerStop(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1062 // AliSysInfo::AddStamp(Form("LRec%s_%d",fgkDetectorName[iDet],eventNr));
1063
1064// // unload calibration data
1065// man->UnloadFromCache(calibPath);
1066// //man->ClearCache();
1067// }
1068
1069// man->SetCacheFlag(origCache);
1070
1071// if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1072// AliError(Form("the following detectors were not found: %s",
1073// detStr.Data()));
1074// if (fStopOnError) return kFALSE;
1075// }
1076
1077 eventNr++;
1078 return kTRUE;
1079}
1080
1081//_____________________________________________________________________________
1082Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
1083{
1084// run the local reconstruction
1085 static Int_t eventNr=0;
1086 AliCodeTimerAuto("")
1087
1088 TString detStr = detectors;
1089 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1090 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1091 AliReconstructor* reconstructor = GetReconstructor(iDet);
1092 if (!reconstructor) continue;
1093 AliLoader* loader = fLoader[iDet];
1094 if (!loader) {
1095 AliWarning(Form("No loader is defined for %s!",fgkDetectorName[iDet]));
1096 continue;
1097 }
1098
1099 // conversion of digits
1100 if (fRawReader && reconstructor->HasDigitConversion()) {
1101 AliInfo(Form("converting raw data digits into root objects for %s",
1102 fgkDetectorName[iDet]));
1103 AliCodeTimerAuto(Form("converting raw data digits into root objects for %s",
1104 fgkDetectorName[iDet]));
1105 loader->LoadDigits("update");
1106 loader->CleanDigits();
1107 loader->MakeDigitsContainer();
1108 TTree* digitsTree = loader->TreeD();
1109 reconstructor->ConvertDigits(fRawReader, digitsTree);
1110 loader->WriteDigits("OVERWRITE");
1111 loader->UnloadDigits();
1112 }
1113
1114 // local reconstruction
1115 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1116 AliCodeTimerAuto(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1117 loader->LoadRecPoints("update");
1118 loader->CleanRecPoints();
1119 loader->MakeRecPointsContainer();
1120 TTree* clustersTree = loader->TreeR();
1121 if (fRawReader && !reconstructor->HasDigitConversion()) {
1122 reconstructor->Reconstruct(fRawReader, clustersTree);
1123 } else {
1124 loader->LoadDigits("read");
1125 TTree* digitsTree = loader->TreeD();
1126 if (!digitsTree) {
1127 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
1128 if (fStopOnError) return kFALSE;
1129 } else {
1130 reconstructor->Reconstruct(digitsTree, clustersTree);
1131 }
1132 loader->UnloadDigits();
1133 }
1134
1135 // do QA in the event loop if requested
1136 if (fQAInLoop) {
1137 AliQADataMaker * qadm = GetQADataMaker(iDet);
1138 if (qadm) {
1139 AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1140 AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1141 if (qadm->IsCycleDone() ) {
1142 qadm->EndOfCycle(AliQA::kRECPOINTS) ;
1143 qadm->EndOfCycle(AliQA::kESDS) ;
1144 qadm->StartOfCycle(AliQA::kRECPOINTS) ;
1145 qadm->StartOfCycle(AliQA::kESDS, "same") ;
1146 }
1147 qadm->Exec(AliQA::kRECPOINTS, clustersTree) ;
1148 AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1149 }
1150 }
1151
1152 loader->WriteRecPoints("OVERWRITE");
1153 loader->UnloadRecPoints();
1154 AliSysInfo::AddStamp(Form("LRec%s_%d",fgkDetectorName[iDet],eventNr), iDet,1,eventNr);
1155 }
1156
1157 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1158 AliError(Form("the following detectors were not found: %s",
1159 detStr.Data()));
1160 if (fStopOnError) return kFALSE;
1161 }
1162 eventNr++;
1163 return kTRUE;
1164}
1165
1166//_____________________________________________________________________________
1167Bool_t AliReconstruction::RunVertexFinder(AliESDEvent*& esd)
1168{
1169// run the barrel tracking
1170
1171 AliCodeTimerAuto("")
1172
1173 AliESDVertex* vertex = NULL;
1174 Double_t vtxPos[3] = {0, 0, 0};
1175 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
1176 TArrayF mcVertex(3);
1177 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
1178 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
1179 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
1180 }
1181
1182 if (fVertexer) {
1183 if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
1184 AliInfo("running the ITS vertex finder");
1185 if (fLoader[0]) fLoader[0]->LoadRecPoints();
1186 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
1187 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
1188 if(!vertex){
1189 AliWarning("Vertex not found");
1190 vertex = new AliESDVertex();
1191 vertex->SetName("default");
1192 }
1193 else {
1194 vertex->SetName("reconstructed");
1195 }
1196
1197 } else {
1198 AliInfo("getting the primary vertex from MC");
1199 vertex = new AliESDVertex(vtxPos, vtxErr);
1200 }
1201
1202 if (vertex) {
1203 vertex->GetXYZ(vtxPos);
1204 vertex->GetSigmaXYZ(vtxErr);
1205 } else {
1206 AliWarning("no vertex reconstructed");
1207 vertex = new AliESDVertex(vtxPos, vtxErr);
1208 }
1209 esd->SetVertex(vertex);
1210 // if SPD multiplicity has been determined, it is stored in the ESD
1211 AliMultiplicity *mult = fVertexer->GetMultiplicity();
1212 if(mult)esd->SetMultiplicity(mult);
1213
1214 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1215 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1216 }
1217 delete vertex;
1218
1219 return kTRUE;
1220}
1221
1222//_____________________________________________________________________________
1223Bool_t AliReconstruction::RunHLTTracking(AliESDEvent*& esd)
1224{
1225// run the HLT barrel tracking
1226
1227 AliCodeTimerAuto("")
1228
1229 if (!fRunLoader) {
1230 AliError("Missing runLoader!");
1231 return kFALSE;
1232 }
1233
1234 AliInfo("running HLT tracking");
1235
1236 // Get a pointer to the HLT reconstructor
1237 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1238 if (!reconstructor) return kFALSE;
1239
1240 // TPC + ITS
1241 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1242 TString detName = fgkDetectorName[iDet];
1243 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1244 reconstructor->SetOption(detName.Data());
1245 AliTracker *tracker = reconstructor->CreateTracker();
1246 if (!tracker) {
1247 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1248 if (fStopOnError) return kFALSE;
1249 continue;
1250 }
1251 Double_t vtxPos[3];
1252 Double_t vtxErr[3]={0.005,0.005,0.010};
1253 const AliESDVertex *vertex = esd->GetVertex();
1254 vertex->GetXYZ(vtxPos);
1255 tracker->SetVertex(vtxPos,vtxErr);
1256 if(iDet != 1) {
1257 fLoader[iDet]->LoadRecPoints("read");
1258 TTree* tree = fLoader[iDet]->TreeR();
1259 if (!tree) {
1260 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1261 return kFALSE;
1262 }
1263 tracker->LoadClusters(tree);
1264 }
1265 if (tracker->Clusters2Tracks(esd) != 0) {
1266 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1267 return kFALSE;
1268 }
1269 if(iDet != 1) {
1270 tracker->UnloadClusters();
1271 }
1272 delete tracker;
1273 }
1274
1275 return kTRUE;
1276}
1277
1278//_____________________________________________________________________________
1279Bool_t AliReconstruction::RunMuonTracking(AliESDEvent*& esd)
1280{
1281// run the muon spectrometer tracking
1282
1283 AliCodeTimerAuto("")
1284
1285 if (!fRunLoader) {
1286 AliError("Missing runLoader!");
1287 return kFALSE;
1288 }
1289 Int_t iDet = 7; // for MUON
1290
1291 AliInfo("is running...");
1292
1293 // Get a pointer to the MUON reconstructor
1294 AliReconstructor *reconstructor = GetReconstructor(iDet);
1295 if (!reconstructor) return kFALSE;
1296
1297
1298 TString detName = fgkDetectorName[iDet];
1299 AliDebug(1, Form("%s tracking", detName.Data()));
1300 AliTracker *tracker = reconstructor->CreateTracker();
1301 if (!tracker) {
1302 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1303 return kFALSE;
1304 }
1305
1306 // read RecPoints
1307 fLoader[iDet]->LoadRecPoints("read");
1308
1309 tracker->LoadClusters(fLoader[iDet]->TreeR());
1310
1311 Int_t rv = tracker->Clusters2Tracks(esd);
1312
1313 if ( rv )
1314 {
1315 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1316 return kFALSE;
1317 }
1318
1319 fLoader[iDet]->UnloadRecPoints();
1320
1321 tracker->UnloadClusters();
1322
1323 delete tracker;
1324
1325 return kTRUE;
1326}
1327
1328
1329//_____________________________________________________________________________
1330Bool_t AliReconstruction::RunTracking(AliESDEvent*& esd)
1331{
1332// run the barrel tracking
1333 static Int_t eventNr=0;
1334 AliCodeTimerAuto("")
1335
1336 AliInfo("running tracking");
1337
1338 //Fill the ESD with the T0 info (will be used by the TOF)
1339 if (fReconstructor[11] && fLoader[11]) {
1340 fLoader[11]->LoadRecPoints("READ");
1341 TTree *treeR = fLoader[11]->TreeR();
1342 GetReconstructor(11)->FillESD((TTree *)NULL,treeR,esd);
1343 }
1344
1345 // pass 1: TPC + ITS inwards
1346 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1347 if (!fTracker[iDet]) continue;
1348 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1349
1350 // load clusters
1351 fLoader[iDet]->LoadRecPoints("read");
1352 AliSysInfo::AddStamp(Form("RLoadCluster%s_%d",fgkDetectorName[iDet],eventNr),iDet,1, eventNr);
1353 TTree* tree = fLoader[iDet]->TreeR();
1354 if (!tree) {
1355 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1356 return kFALSE;
1357 }
1358 fTracker[iDet]->LoadClusters(tree);
1359 AliSysInfo::AddStamp(Form("TLoadCluster%s_%d",fgkDetectorName[iDet],eventNr), iDet,2, eventNr);
1360 // run tracking
1361 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1362 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1363 return kFALSE;
1364 }
1365 if (fCheckPointLevel > 1) {
1366 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1367 }
1368 // preliminary PID in TPC needed by the ITS tracker
1369 if (iDet == 1) {
1370 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1371 AliESDpid::MakePID(esd);
1372 }
1373 AliSysInfo::AddStamp(Form("Tracking0%s_%d",fgkDetectorName[iDet],eventNr), iDet,3,eventNr);
1374 }
1375
1376 // pass 2: ALL backwards
1377 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1378 if (!fTracker[iDet]) continue;
1379 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1380
1381 // load clusters
1382 if (iDet > 1) { // all except ITS, TPC
1383 TTree* tree = NULL;
1384 fLoader[iDet]->LoadRecPoints("read");
1385 AliSysInfo::AddStamp(Form("RLoadCluster0%s_%d",fgkDetectorName[iDet],eventNr), iDet,1, eventNr);
1386 tree = fLoader[iDet]->TreeR();
1387 if (!tree) {
1388 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1389 return kFALSE;
1390 }
1391 fTracker[iDet]->LoadClusters(tree);
1392 AliSysInfo::AddStamp(Form("TLoadCluster0%s_%d",fgkDetectorName[iDet],eventNr), iDet,2, eventNr);
1393 }
1394
1395 // run tracking
1396 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1397 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1398 // return kFALSE;
1399 }
1400 if (fCheckPointLevel > 1) {
1401 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1402 }
1403
1404 // unload clusters
1405 if (iDet > 2) { // all except ITS, TPC, TRD
1406 fTracker[iDet]->UnloadClusters();
1407 fLoader[iDet]->UnloadRecPoints();
1408 }
1409 // updated PID in TPC needed by the ITS tracker -MI
1410 if (iDet == 1) {
1411 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1412 AliESDpid::MakePID(esd);
1413 }
1414 AliSysInfo::AddStamp(Form("Tracking1%s_%d",fgkDetectorName[iDet],eventNr), iDet,3, eventNr);
1415 }
1416
1417 // write space-points to the ESD in case alignment data output
1418 // is switched on
1419 if (fWriteAlignmentData)
1420 WriteAlignmentData(esd);
1421
1422 // pass 3: TRD + TPC + ITS refit inwards
1423 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1424 if (!fTracker[iDet]) continue;
1425 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1426
1427 // run tracking
1428 if (fTracker[iDet]->RefitInward(esd) != 0) {
1429 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1430 // return kFALSE;
1431 }
1432 if (fCheckPointLevel > 1) {
1433 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1434 }
1435 AliSysInfo::AddStamp(Form("Tracking2%s_%d",fgkDetectorName[iDet],eventNr), iDet,3, eventNr);
1436 // unload clusters
1437 fTracker[iDet]->UnloadClusters();
1438 AliSysInfo::AddStamp(Form("TUnloadCluster%s_%d",fgkDetectorName[iDet],eventNr), iDet,4, eventNr);
1439 fLoader[iDet]->UnloadRecPoints();
1440 AliSysInfo::AddStamp(Form("RUnloadCluster%s_%d",fgkDetectorName[iDet],eventNr), iDet,5, eventNr);
1441 }
1442 //
1443 // Propagate track to the vertex - if not done by ITS
1444 //
1445 Int_t ntracks = esd->GetNumberOfTracks();
1446 for (Int_t itrack=0; itrack<ntracks; itrack++){
1447 const Double_t kRadius = 3; // beam pipe radius
1448 const Double_t kMaxStep = 5; // max step
1449 const Double_t kMaxD = 123456; // max distance to prim vertex
1450 Double_t fieldZ = AliTracker::GetBz(); //
1451 AliESDtrack * track = esd->GetTrack(itrack);
1452 if (!track) continue;
1453 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1454 AliTracker::PropagateTrackTo(track,kRadius,track->GetMass(),kMaxStep,kTRUE);
1455 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1456 }
1457 eventNr++;
1458 return kTRUE;
1459}
1460
1461//_____________________________________________________________________________
1462Bool_t AliReconstruction::CleanESD(AliESDEvent *esd){
1463 //
1464 // Remove the data which are not needed for the physics analysis.
1465 //
1466
1467 Int_t nTracks=esd->GetNumberOfTracks();
1468 Int_t nV0s=esd->GetNumberOfV0s();
1469 AliInfo
1470 (Form("Number of ESD tracks and V0s before cleaning: %d %d",nTracks,nV0s));
1471
1472 Float_t cleanPars[]={fV0DCAmax,fV0CsPmin,fDmax,fZmax};
1473 Bool_t rc=esd->Clean(cleanPars);
1474
1475 nTracks=esd->GetNumberOfTracks();
1476 nV0s=esd->GetNumberOfV0s();
1477 AliInfo
1478 (Form("Number of ESD tracks and V0s after cleaning %d",nTracks,nV0s));
1479
1480 return rc;
1481}
1482
1483//_____________________________________________________________________________
1484Bool_t AliReconstruction::FillESD(AliESDEvent*& esd, const TString& detectors)
1485{
1486// fill the event summary data
1487
1488 AliCodeTimerAuto("")
1489 static Int_t eventNr=0;
1490 TString detStr = detectors;
1491
1492 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1493 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1494 AliReconstructor* reconstructor = GetReconstructor(iDet);
1495 if (!reconstructor) continue;
1496 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1497 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1498 TTree* clustersTree = NULL;
1499 if (fLoader[iDet]) {
1500 fLoader[iDet]->LoadRecPoints("read");
1501 clustersTree = fLoader[iDet]->TreeR();
1502 if (!clustersTree) {
1503 AliError(Form("Can't get the %s clusters tree",
1504 fgkDetectorName[iDet]));
1505 if (fStopOnError) return kFALSE;
1506 }
1507 }
1508 if (fRawReader && !reconstructor->HasDigitConversion()) {
1509 reconstructor->FillESD(fRawReader, clustersTree, esd);
1510 } else {
1511 TTree* digitsTree = NULL;
1512 if (fLoader[iDet]) {
1513 fLoader[iDet]->LoadDigits("read");
1514 digitsTree = fLoader[iDet]->TreeD();
1515 if (!digitsTree) {
1516 AliError(Form("Can't get the %s digits tree",
1517 fgkDetectorName[iDet]));
1518 if (fStopOnError) return kFALSE;
1519 }
1520 }
1521 reconstructor->FillESD(digitsTree, clustersTree, esd);
1522 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1523 }
1524 if (fLoader[iDet]) {
1525 fLoader[iDet]->UnloadRecPoints();
1526 }
1527
1528 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1529 }
1530 }
1531
1532 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1533 AliError(Form("the following detectors were not found: %s",
1534 detStr.Data()));
1535 if (fStopOnError) return kFALSE;
1536 }
1537 AliSysInfo::AddStamp(Form("FillESD%d",eventNr), 0,1, eventNr);
1538 eventNr++;
1539 return kTRUE;
1540}
1541
1542//_____________________________________________________________________________
1543Bool_t AliReconstruction::FillTriggerESD(AliESDEvent*& esd)
1544{
1545 // Reads the trigger decision which is
1546 // stored in Trigger.root file and fills
1547 // the corresponding esd entries
1548
1549 AliCodeTimerAuto("")
1550
1551 AliInfo("Filling trigger information into the ESD");
1552
1553 if (fRawReader) {
1554 AliCTPRawStream input(fRawReader);
1555 if (!input.Next()) {
1556 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1557 return kFALSE;
1558 }
1559 esd->SetTriggerMask(input.GetClassMask());
1560 esd->SetTriggerCluster(input.GetClusterMask());
1561 }
1562 else {
1563 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1564 if (runloader) {
1565 if (!runloader->LoadTrigger()) {
1566 AliCentralTrigger *aCTP = runloader->GetTrigger();
1567 esd->SetTriggerMask(aCTP->GetClassMask());
1568 esd->SetTriggerCluster(aCTP->GetClusterMask());
1569 }
1570 else {
1571 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1572 return kFALSE;
1573 }
1574 }
1575 else {
1576 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1577 return kFALSE;
1578 }
1579 }
1580
1581 return kTRUE;
1582}
1583
1584
1585
1586
1587
1588//_____________________________________________________________________________
1589Bool_t AliReconstruction::FillRawEventHeaderESD(AliESDEvent*& esd)
1590{
1591 //
1592 // Filling information from RawReader Header
1593 //
1594
1595 AliInfo("Filling information from RawReader Header");
1596 esd->SetBunchCrossNumber(0);
1597 esd->SetOrbitNumber(0);
1598 esd->SetPeriodNumber(0);
1599 esd->SetTimeStamp(0);
1600 esd->SetEventType(0);
1601 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1602 if (eventHeader){
1603
1604 const UInt_t *id = eventHeader->GetP("Id");
1605 esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1606 esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1607 esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1608
1609 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1610 esd->SetEventType((eventHeader->Get("Type")));
1611 }
1612
1613 return kTRUE;
1614}
1615
1616
1617//_____________________________________________________________________________
1618Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1619{
1620// check whether detName is contained in detectors
1621// if yes, it is removed from detectors
1622
1623 // check if all detectors are selected
1624 if ((detectors.CompareTo("ALL") == 0) ||
1625 detectors.BeginsWith("ALL ") ||
1626 detectors.EndsWith(" ALL") ||
1627 detectors.Contains(" ALL ")) {
1628 detectors = "ALL";
1629 return kTRUE;
1630 }
1631
1632 // search for the given detector
1633 Bool_t result = kFALSE;
1634 if ((detectors.CompareTo(detName) == 0) ||
1635 detectors.BeginsWith(detName+" ") ||
1636 detectors.EndsWith(" "+detName) ||
1637 detectors.Contains(" "+detName+" ")) {
1638 detectors.ReplaceAll(detName, "");
1639 result = kTRUE;
1640 }
1641
1642 // clean up the detectors string
1643 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1644 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1645 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1646
1647 return result;
1648}
1649
1650//_____________________________________________________________________________
1651Bool_t AliReconstruction::InitRunLoader()
1652{
1653// get or create the run loader
1654
1655 if (gAlice) delete gAlice;
1656 gAlice = NULL;
1657
1658 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1659 // load all base libraries to get the loader classes
1660 TString libs = gSystem->GetLibraries();
1661 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1662 TString detName = fgkDetectorName[iDet];
1663 if (detName == "HLT") continue;
1664 if (libs.Contains("lib" + detName + "base.so")) continue;
1665 gSystem->Load("lib" + detName + "base.so");
1666 }
1667 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1668 if (!fRunLoader) {
1669 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1670 CleanUp();
1671 return kFALSE;
1672 }
1673 fRunLoader->CdGAFile();
1674 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1675 if (fRunLoader->LoadgAlice() == 0) {
1676 gAlice = fRunLoader->GetAliRun();
1677 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1678 }
1679 }
1680 if (!gAlice && !fRawReader) {
1681 AliError(Form("no gAlice object found in file %s",
1682 fGAliceFileName.Data()));
1683 CleanUp();
1684 return kFALSE;
1685 }
1686
1687 //PH This is a temporary fix to give access to the kinematics
1688 //PH that is needed for the labels of ITS clusters
1689 fRunLoader->LoadHeader();
1690 fRunLoader->LoadKinematics();
1691
1692 } else { // galice.root does not exist
1693 if (!fRawReader) {
1694 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1695 CleanUp();
1696 return kFALSE;
1697 }
1698 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1699 AliConfig::GetDefaultEventFolderName(),
1700 "recreate");
1701 if (!fRunLoader) {
1702 AliError(Form("could not create run loader in file %s",
1703 fGAliceFileName.Data()));
1704 CleanUp();
1705 return kFALSE;
1706 }
1707 fRunLoader->MakeTree("E");
1708 Int_t iEvent = 0;
1709 while (fRawReader->NextEvent()) {
1710 fRunLoader->SetEventNumber(iEvent);
1711 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1712 iEvent, iEvent);
1713 fRunLoader->MakeTree("H");
1714 fRunLoader->TreeE()->Fill();
1715 iEvent++;
1716 }
1717 fRawReader->RewindEvents();
1718 if (fNumberOfEventsPerFile > 0)
1719 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1720 else
1721 fRunLoader->SetNumberOfEventsPerFile(iEvent);
1722 fRunLoader->WriteHeader("OVERWRITE");
1723 fRunLoader->CdGAFile();
1724 fRunLoader->Write(0, TObject::kOverwrite);
1725// AliTracker::SetFieldMap(???);
1726 }
1727
1728 return kTRUE;
1729}
1730
1731//_____________________________________________________________________________
1732AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1733{
1734// get the reconstructor object and the loader for a detector
1735
1736 if (fReconstructor[iDet]) return fReconstructor[iDet];
1737
1738 // load the reconstructor object
1739 TPluginManager* pluginManager = gROOT->GetPluginManager();
1740 TString detName = fgkDetectorName[iDet];
1741 TString recName = "Ali" + detName + "Reconstructor";
1742 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1743
1744 AliReconstructor* reconstructor = NULL;
1745 // first check if a plugin is defined for the reconstructor
1746 TPluginHandler* pluginHandler =
1747 pluginManager->FindHandler("AliReconstructor", detName);
1748 // if not, add a plugin for it
1749 if (!pluginHandler) {
1750 AliDebug(1, Form("defining plugin for %s", recName.Data()));
1751 TString libs = gSystem->GetLibraries();
1752 if (libs.Contains("lib" + detName + "base.so") ||
1753 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1754 pluginManager->AddHandler("AliReconstructor", detName,
1755 recName, detName + "rec", recName + "()");
1756 } else {
1757 pluginManager->AddHandler("AliReconstructor", detName,
1758 recName, detName, recName + "()");
1759 }
1760 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1761 }
1762 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1763 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1764 }
1765 if (reconstructor) {
1766 TObject* obj = fOptions.FindObject(detName.Data());
1767 if (obj) reconstructor->SetOption(obj->GetTitle());
1768 reconstructor->Init();
1769 fReconstructor[iDet] = reconstructor;
1770 }
1771
1772 // get or create the loader
1773 if (detName != "HLT") {
1774 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1775 if (!fLoader[iDet]) {
1776 AliConfig::Instance()
1777 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1778 detName, detName);
1779 // first check if a plugin is defined for the loader
1780 pluginHandler =
1781 pluginManager->FindHandler("AliLoader", detName);
1782 // if not, add a plugin for it
1783 if (!pluginHandler) {
1784 TString loaderName = "Ali" + detName + "Loader";
1785 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1786 pluginManager->AddHandler("AliLoader", detName,
1787 loaderName, detName + "base",
1788 loaderName + "(const char*, TFolder*)");
1789 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1790 }
1791 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1792 fLoader[iDet] =
1793 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1794 fRunLoader->GetEventFolder());
1795 }
1796 if (!fLoader[iDet]) { // use default loader
1797 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1798 }
1799 if (!fLoader[iDet]) {
1800 AliWarning(Form("couldn't get loader for %s", detName.Data()));
1801 if (fStopOnError) return NULL;
1802 } else {
1803 fRunLoader->AddLoader(fLoader[iDet]);
1804 fRunLoader->CdGAFile();
1805 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1806 fRunLoader->Write(0, TObject::kOverwrite);
1807 }
1808 }
1809 }
1810
1811 return reconstructor;
1812}
1813
1814//_____________________________________________________________________________
1815Bool_t AliReconstruction::CreateVertexer()
1816{
1817// create the vertexer
1818
1819 fVertexer = NULL;
1820 AliReconstructor* itsReconstructor = GetReconstructor(0);
1821 if (itsReconstructor) {
1822 fVertexer = itsReconstructor->CreateVertexer();
1823 }
1824 if (!fVertexer) {
1825 AliWarning("couldn't create a vertexer for ITS");
1826 if (fStopOnError) return kFALSE;
1827 }
1828
1829 return kTRUE;
1830}
1831
1832//_____________________________________________________________________________
1833Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1834{
1835// create the trackers
1836
1837 TString detStr = detectors;
1838 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1839 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1840 AliReconstructor* reconstructor = GetReconstructor(iDet);
1841 if (!reconstructor) continue;
1842 TString detName = fgkDetectorName[iDet];
1843 if (detName == "HLT") {
1844 fRunHLTTracking = kTRUE;
1845 continue;
1846 }
1847 if (detName == "MUON") {
1848 fRunMuonTracking = kTRUE;
1849 continue;
1850 }
1851
1852
1853 fTracker[iDet] = reconstructor->CreateTracker();
1854 if (!fTracker[iDet] && (iDet < 7)) {
1855 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1856 if (fStopOnError) return kFALSE;
1857 }
1858 AliSysInfo::AddStamp(Form("LTracker%s",fgkDetectorName[iDet]), iDet,0);
1859 }
1860
1861 return kTRUE;
1862}
1863
1864//_____________________________________________________________________________
1865void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1866{
1867// delete trackers and the run loader and close and delete the file
1868
1869 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1870 delete fReconstructor[iDet];
1871 fReconstructor[iDet] = NULL;
1872 fLoader[iDet] = NULL;
1873 delete fTracker[iDet];
1874 fTracker[iDet] = NULL;
1875// delete fQADataMaker[iDet];
1876// fQADataMaker[iDet] = NULL;
1877 }
1878 delete fVertexer;
1879 fVertexer = NULL;
1880
1881 if(!(AliCDBManager::Instance()->GetCacheFlag())) {
1882 delete fDiamondProfile;
1883 fDiamondProfile = NULL;
1884 }
1885
1886 delete fGRPList;
1887 fGRPList = NULL;
1888
1889 delete fRunLoader;
1890 fRunLoader = NULL;
1891 delete fRawReader;
1892 fRawReader = NULL;
1893
1894 if (file) {
1895 file->Close();
1896 delete file;
1897 }
1898
1899 if (fileOld) {
1900 fileOld->Close();
1901 delete fileOld;
1902 gSystem->Unlink("AliESDs.old.root");
1903 }
1904}
1905
1906//_____________________________________________________________________________
1907
1908Bool_t AliReconstruction::ReadESD(AliESDEvent*& esd, const char* recStep) const
1909{
1910// read the ESD event from a file
1911
1912 if (!esd) return kFALSE;
1913 char fileName[256];
1914 sprintf(fileName, "ESD_%d.%d_%s.root",
1915 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1916 if (gSystem->AccessPathName(fileName)) return kFALSE;
1917
1918 AliInfo(Form("reading ESD from file %s", fileName));
1919 AliDebug(1, Form("reading ESD from file %s", fileName));
1920 TFile* file = TFile::Open(fileName);
1921 if (!file || !file->IsOpen()) {
1922 AliError(Form("opening %s failed", fileName));
1923 delete file;
1924 return kFALSE;
1925 }
1926
1927 gROOT->cd();
1928 delete esd;
1929 esd = (AliESDEvent*) file->Get("ESD");
1930 file->Close();
1931 delete file;
1932 return kTRUE;
1933
1934}
1935
1936
1937
1938//_____________________________________________________________________________
1939void AliReconstruction::WriteESD(AliESDEvent* esd, const char* recStep) const
1940{
1941// write the ESD event to a file
1942
1943 if (!esd) return;
1944 char fileName[256];
1945 sprintf(fileName, "ESD_%d.%d_%s.root",
1946 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1947
1948 AliDebug(1, Form("writing ESD to file %s", fileName));
1949 TFile* file = TFile::Open(fileName, "recreate");
1950 if (!file || !file->IsOpen()) {
1951 AliError(Form("opening %s failed", fileName));
1952 } else {
1953 esd->Write("ESD");
1954 file->Close();
1955 }
1956 delete file;
1957}
1958
1959
1960
1961
1962
1963//_____________________________________________________________________________
1964void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
1965{
1966 // write all files from the given esd file to an aod file
1967
1968 // create an AliAOD object
1969 AliAODEvent *aod = new AliAODEvent();
1970 aod->CreateStdContent();
1971
1972 // go to the file
1973 aodFile->cd();
1974
1975 // create the tree
1976 TTree *aodTree = new TTree("aodTree", "AliAOD tree");
1977 aodTree->Branch(aod->GetList());
1978
1979 // connect to ESD
1980 TTree *t = (TTree*) esdFile->Get("esdTree");
1981 AliESDEvent *esd = new AliESDEvent();
1982 esd->ReadFromTree(t);
1983
1984 Int_t nEvents = t->GetEntries();
1985
1986 // set arrays and pointers
1987 Float_t posF[3];
1988 Double_t pos[3];
1989 Double_t p[3];
1990 Double_t p_pos[3];
1991 Double_t p_neg[3];
1992 Double_t covVtx[6];
1993 Double_t covTr[21];
1994 Double_t pid[10];
1995
1996 // loop over events and fill them
1997 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
1998 //cout << "event: " << iEvent << endl;
1999 t->GetEntry(iEvent);
2000
2001 // Multiplicity information needed by the header (to be revised!)
2002 Int_t nTracks = esd->GetNumberOfTracks();
2003 Int_t nPosTracks = 0;
2004 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2005 if (esd->GetTrack(iTrack)->Charge()> 0) nPosTracks++;
2006
2007 // Access the header
2008 AliAODHeader *header = aod->GetHeader();
2009
2010 // fill the header
2011 header->SetRunNumber (esd->GetRunNumber() );
2012 header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
2013 header->SetOrbitNumber (esd->GetOrbitNumber() );
2014 header->SetPeriodNumber (esd->GetPeriodNumber() );
2015 header->SetTriggerMask (esd->GetTriggerMask() );
2016 header->SetTriggerCluster (esd->GetTriggerCluster() );
2017 header->SetEventType (esd->GetEventType() );
2018 header->SetMagneticField (esd->GetMagneticField() );
2019 header->SetZDCN1Energy (esd->GetZDCN1Energy() );
2020 header->SetZDCP1Energy (esd->GetZDCP1Energy() );
2021 header->SetZDCN2Energy (esd->GetZDCN2Energy() );
2022 header->SetZDCP2Energy (esd->GetZDCP2Energy() );
2023 header->SetZDCEMEnergy (esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
2024 header->SetRefMultiplicity (nTracks);
2025 header->SetRefMultiplicityPos(nPosTracks);
2026 header->SetRefMultiplicityNeg(nTracks - nPosTracks);
2027 header->SetMuonMagFieldScale(-999.); // FIXME
2028 header->SetCentrality(-999.); // FIXME
2029
2030 Int_t nV0s = esd->GetNumberOfV0s();
2031 Int_t nCascades = esd->GetNumberOfCascades();
2032 Int_t nKinks = esd->GetNumberOfKinks();
2033 Int_t nVertices = nV0s + 2*nCascades /*could lead to two vertices, one V0 and the Xi */+ nKinks + 1 /* = prim. vtx*/;
2034 Int_t nJets = 0;
2035 Int_t nCaloClus = esd->GetNumberOfCaloClusters();
2036 Int_t nFmdClus = 0;
2037 Int_t nPmdClus = esd->GetNumberOfPmdTracks();
2038
2039 aod->ResetStd(nTracks, nVertices, nV0s+nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
2040
2041 // Array to take into account the tracks already added to the AOD
2042 Bool_t * usedTrack = NULL;
2043 if (nTracks>0) {
2044 usedTrack = new Bool_t[nTracks];
2045 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2046 }
2047 // Array to take into account the V0s already added to the AOD
2048 Bool_t * usedV0 = NULL;
2049 if (nV0s>0) {
2050 usedV0 = new Bool_t[nV0s];
2051 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2052 }
2053 // Array to take into account the kinks already added to the AOD
2054 Bool_t * usedKink = NULL;
2055 if (nKinks>0) {
2056 usedKink = new Bool_t[nKinks];
2057 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2058 }
2059
2060 // Access to the AOD container of vertices
2061 TClonesArray &vertices = *(aod->GetVertices());
2062 Int_t jVertices=0;
2063
2064 // Access to the AOD container of tracks
2065 TClonesArray &tracks = *(aod->GetTracks());
2066 Int_t jTracks=0;
2067
2068 // Access to the AOD container of V0s
2069 TClonesArray &V0s = *(aod->GetV0s());
2070 Int_t jV0s=0;
2071
2072 // Add primary vertex. The primary tracks will be defined
2073 // after the loops on the composite objects (V0, cascades, kinks)
2074 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2075
2076 vtx->GetXYZ(pos); // position
2077 vtx->GetCovMatrix(covVtx); //covariance matrix
2078
2079 AliAODVertex * primary = new(vertices[jVertices++])
2080 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
2081
2082
2083 AliAODTrack *aodTrack = 0x0;
2084
2085 // Create vertices starting from the most complex objects
2086
2087 // Cascades
2088 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2089 AliESDcascade *cascade = esd->GetCascade(nCascade);
2090
2091 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2092 cascade->GetPosCovXi(covVtx);
2093
2094 // Add the cascade vertex
2095 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2096 covVtx,
2097 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2098 primary,
2099 nCascade,
2100 AliAODVertex::kCascade);
2101
2102 primary->AddDaughter(vcascade); // the cascade 'particle' (represented by a vertex) is added as a daughter to the primary vertex
2103
2104 // Add the V0 from the cascade. The ESD class have to be optimized...
2105 // Now we have to search for the corresponding V0 in the list of V0s
2106 // using the indeces of the positive and negative tracks
2107
2108 Int_t posFromV0 = cascade->GetPindex();
2109 Int_t negFromV0 = cascade->GetNindex();
2110
2111
2112 AliESDv0 * v0 = 0x0;
2113 Int_t indV0 = -1;
2114
2115 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2116
2117 v0 = esd->GetV0(iV0);
2118 Int_t posV0 = v0->GetPindex();
2119 Int_t negV0 = v0->GetNindex();
2120
2121 if (posV0==posFromV0 && negV0==negFromV0) {
2122 indV0 = iV0;
2123 break;
2124 }
2125 }
2126
2127 AliAODVertex * vV0FromCascade = 0x0;
2128
2129 if (indV0>-1 && !usedV0[indV0]) {
2130
2131 // the V0 exists in the array of V0s and is not used
2132
2133 usedV0[indV0] = kTRUE;
2134
2135 v0->GetXYZ(pos[0], pos[1], pos[2]);
2136 v0->GetPosCov(covVtx);
2137
2138 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2139 covVtx,
2140 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2141 vcascade,
2142 indV0,
2143 AliAODVertex::kV0);
2144 } else {
2145
2146 // the V0 doesn't exist in the array of V0s or was used
2147 cerr << "Error: event " << iEvent << " cascade " << nCascade
2148 << " The V0 " << indV0
2149 << " doesn't exist in the array of V0s or was used!" << endl;
2150
2151 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2152 cascade->GetPosCov(covVtx);
2153
2154 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2155 covVtx,
2156 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2157 vcascade,
2158 indV0,
2159 AliAODVertex::kV0);
2160 vcascade->AddDaughter(vV0FromCascade);
2161
2162 }
2163
2164 // Add the positive tracks from the V0
2165
2166 if (! usedTrack[posFromV0]) {
2167
2168 usedTrack[posFromV0] = kTRUE;
2169
2170 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2171 esdTrack->GetPxPyPz(p_pos);
2172 esdTrack->GetXYZ(pos);
2173 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2174 esdTrack->GetESDpid(pid);
2175
2176 vV0FromCascade->AddDaughter(aodTrack =
2177 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2178 esdTrack->GetLabel(),
2179 p_pos,
2180 kTRUE,
2181 pos,
2182 kFALSE,
2183 covTr,
2184 (Short_t)esdTrack->Charge(),
2185 esdTrack->GetITSClusterMap(),
2186 pid,
2187 vV0FromCascade,
2188 kTRUE, // check if this is right
2189 kFALSE, // check if this is right
2190 AliAODTrack::kSecondary)
2191 );
2192 aodTrack->ConvertAliPIDtoAODPID();
2193 }
2194 else {
2195 cerr << "Error: event " << iEvent << " cascade " << nCascade
2196 << " track " << posFromV0 << " has already been used!" << endl;
2197 }
2198
2199 // Add the negative tracks from the V0
2200
2201 if (!usedTrack[negFromV0]) {
2202
2203 usedTrack[negFromV0] = kTRUE;
2204
2205 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2206 esdTrack->GetPxPyPz(p_neg);
2207 esdTrack->GetXYZ(pos);
2208 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2209 esdTrack->GetESDpid(pid);
2210
2211 vV0FromCascade->AddDaughter(aodTrack =
2212 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2213 esdTrack->GetLabel(),
2214 p_neg,
2215 kTRUE,
2216 pos,
2217 kFALSE,
2218 covTr,
2219 (Short_t)esdTrack->Charge(),
2220 esdTrack->GetITSClusterMap(),
2221 pid,
2222 vV0FromCascade,
2223 kTRUE, // check if this is right
2224 kFALSE, // check if this is right
2225 AliAODTrack::kSecondary)
2226 );
2227 aodTrack->ConvertAliPIDtoAODPID();
2228 }
2229 else {
2230 cerr << "Error: event " << iEvent << " cascade " << nCascade
2231 << " track " << negFromV0 << " has already been used!" << endl;
2232 }
2233
2234 // add it to the V0 array as well
2235 Double_t d0[2] = { -999., -99.};
2236 // counting is probably wrong
2237 new(V0s[jV0s++]) AliAODv0(vV0FromCascade, -999., -99., p_pos, p_neg, d0); // to be refined
2238
2239 // Add the bachelor track from the cascade
2240
2241 Int_t bachelor = cascade->GetBindex();
2242
2243 if(!usedTrack[bachelor]) {
2244
2245 usedTrack[bachelor] = kTRUE;
2246
2247 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2248 esdTrack->GetPxPyPz(p);
2249 esdTrack->GetXYZ(pos);
2250 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2251 esdTrack->GetESDpid(pid);
2252
2253 vcascade->AddDaughter(aodTrack =
2254 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2255 esdTrack->GetLabel(),
2256 p,
2257 kTRUE,
2258 pos,
2259 kFALSE,
2260 covTr,
2261 (Short_t)esdTrack->Charge(),
2262 esdTrack->GetITSClusterMap(),
2263 pid,
2264 vcascade,
2265 kTRUE, // check if this is right
2266 kFALSE, // check if this is right
2267 AliAODTrack::kSecondary)
2268 );
2269 aodTrack->ConvertAliPIDtoAODPID();
2270 }
2271 else {
2272 cerr << "Error: event " << iEvent << " cascade " << nCascade
2273 << " track " << bachelor << " has already been used!" << endl;
2274 }
2275
2276 // Add the primary track of the cascade (if any)
2277
2278 } // end of the loop on cascades
2279
2280 // V0s
2281
2282 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2283
2284 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2285
2286 AliESDv0 *v0 = esd->GetV0(nV0);
2287
2288 v0->GetXYZ(pos[0], pos[1], pos[2]);
2289 v0->GetPosCov(covVtx);
2290
2291 AliAODVertex * vV0 =
2292 new(vertices[jVertices++]) AliAODVertex(pos,
2293 covVtx,
2294 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2295 primary,
2296 nV0,
2297 AliAODVertex::kV0);
2298 primary->AddDaughter(vV0);
2299
2300 Int_t posFromV0 = v0->GetPindex();
2301 Int_t negFromV0 = v0->GetNindex();
2302
2303 // Add the positive tracks from the V0
2304
2305 if (!usedTrack[posFromV0]) {
2306
2307 usedTrack[posFromV0] = kTRUE;
2308
2309 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2310 esdTrack->GetPxPyPz(p_pos);
2311 esdTrack->GetXYZ(pos);
2312 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2313 esdTrack->GetESDpid(pid);
2314
2315 vV0->AddDaughter(aodTrack =
2316 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2317 esdTrack->GetLabel(),
2318 p_pos,
2319 kTRUE,
2320 pos,
2321 kFALSE,
2322 covTr,
2323 (Short_t)esdTrack->Charge(),
2324 esdTrack->GetITSClusterMap(),
2325 pid,
2326 vV0,
2327 kTRUE, // check if this is right
2328 kFALSE, // check if this is right
2329 AliAODTrack::kSecondary)
2330 );
2331 aodTrack->ConvertAliPIDtoAODPID();
2332 }
2333 else {
2334 cerr << "Error: event " << iEvent << " V0 " << nV0
2335 << " track " << posFromV0 << " has already been used!" << endl;
2336 }
2337
2338 // Add the negative tracks from the V0
2339
2340 if (!usedTrack[negFromV0]) {
2341
2342 usedTrack[negFromV0] = kTRUE;
2343
2344 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2345 esdTrack->GetPxPyPz(p_neg);
2346 esdTrack->GetXYZ(pos);
2347 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2348 esdTrack->GetESDpid(pid);
2349
2350 vV0->AddDaughter(aodTrack =
2351 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2352 esdTrack->GetLabel(),
2353 p_neg,
2354 kTRUE,
2355 pos,
2356 kFALSE,
2357 covTr,
2358 (Short_t)esdTrack->Charge(),
2359 esdTrack->GetITSClusterMap(),
2360 pid,
2361 vV0,
2362 kTRUE, // check if this is right
2363 kFALSE, // check if this is right
2364 AliAODTrack::kSecondary)
2365 );
2366 aodTrack->ConvertAliPIDtoAODPID();
2367 }
2368 else {
2369 cerr << "Error: event " << iEvent << " V0 " << nV0
2370 << " track " << negFromV0 << " has already been used!" << endl;
2371 }
2372
2373 // add it to the V0 array as well
2374 Double_t d0[2] = { 999., 99.};
2375 new(V0s[jV0s++]) AliAODv0(vV0, 999., 99., p_pos, p_neg, d0); // to be refined
2376 }
2377 V0s.Expand(jV0s);
2378 // end of the loop on V0s
2379
2380 // Kinks: it is a big mess the access to the information in the kinks
2381 // The loop is on the tracks in order to find the mother and daugther of each kink
2382
2383
2384 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2385
2386 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2387
2388 Int_t ikink = esdTrack->GetKinkIndex(0);
2389
2390 if (ikink) {
2391 // Negative kink index: mother, positive: daughter
2392
2393 // Search for the second track of the kink
2394
2395 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2396
2397 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2398
2399 Int_t jkink = esdTrack1->GetKinkIndex(0);
2400
2401 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2402
2403 // The two tracks are from the same kink
2404
2405 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2406
2407 Int_t imother = -1;
2408 Int_t idaughter = -1;
2409
2410 if (ikink<0 && jkink>0) {
2411
2412 imother = iTrack;
2413 idaughter = jTrack;
2414 }
2415 else if (ikink>0 && jkink<0) {
2416
2417 imother = jTrack;
2418 idaughter = iTrack;
2419 }
2420 else {
2421 cerr << "Error: Wrong combination of kink indexes: "
2422 << ikink << " " << jkink << endl;
2423 continue;
2424 }
2425
2426 // Add the mother track
2427
2428 AliAODTrack * mother = NULL;
2429
2430 if (!usedTrack[imother]) {
2431
2432 usedTrack[imother] = kTRUE;
2433
2434 AliESDtrack *esdTrack = esd->GetTrack(imother);
2435 esdTrack->GetPxPyPz(p);
2436 esdTrack->GetXYZ(pos);
2437 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2438 esdTrack->GetESDpid(pid);
2439
2440 mother =
2441 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2442 esdTrack->GetLabel(),
2443 p,
2444 kTRUE,
2445 pos,
2446 kFALSE,
2447 covTr,
2448 (Short_t)esdTrack->Charge(),
2449 esdTrack->GetITSClusterMap(),
2450 pid,
2451 primary,
2452 kTRUE, // check if this is right
2453 kTRUE, // check if this is right
2454 AliAODTrack::kPrimary);
2455 primary->AddDaughter(mother);
2456 mother->ConvertAliPIDtoAODPID();
2457 }
2458 else {
2459 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2460 << " track " << imother << " has already been used!" << endl;
2461 }
2462
2463 // Add the kink vertex
2464 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2465
2466 AliAODVertex * vkink =
2467 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2468 NULL,
2469 0.,
2470 mother,
2471 esdTrack->GetID(), // This is the track ID of the mother's track!
2472 AliAODVertex::kKink);
2473 // Add the daughter track
2474
2475 AliAODTrack * daughter = NULL;
2476
2477 if (!usedTrack[idaughter]) {
2478
2479 usedTrack[idaughter] = kTRUE;
2480
2481 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2482 esdTrack->GetPxPyPz(p);
2483 esdTrack->GetXYZ(pos);
2484 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2485 esdTrack->GetESDpid(pid);
2486
2487 daughter =
2488 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2489 esdTrack->GetLabel(),
2490 p,
2491 kTRUE,
2492 pos,
2493 kFALSE,
2494 covTr,
2495 (Short_t)esdTrack->Charge(),
2496 esdTrack->GetITSClusterMap(),
2497 pid,
2498 vkink,
2499 kTRUE, // check if this is right
2500 kTRUE, // check if this is right
2501 AliAODTrack::kPrimary);
2502 vkink->AddDaughter(daughter);
2503 daughter->ConvertAliPIDtoAODPID();
2504 }
2505 else {
2506 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2507 << " track " << idaughter << " has already been used!" << endl;
2508 }
2509 }
2510 }
2511 }
2512 }
2513 vertices.Expand(jVertices);
2514
2515 // Tracks (primary and orphan)
2516 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2517
2518 if (usedTrack[nTrack]) continue;
2519
2520 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2521 esdTrack->GetPxPyPz(p);
2522 esdTrack->GetXYZ(pos);
2523 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2524 esdTrack->GetESDpid(pid);
2525
2526 Float_t impactXY, impactZ;
2527
2528 esdTrack->GetImpactParameters(impactXY,impactZ);
2529
2530 if (impactXY<3.) {
2531 // track inside the beam pipe
2532
2533 primary->AddDaughter(aodTrack =
2534 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2535 esdTrack->GetLabel(),
2536 p,
2537 kTRUE,
2538 pos,
2539 kFALSE,
2540 covTr,
2541 (Short_t)esdTrack->Charge(),
2542 esdTrack->GetITSClusterMap(),
2543 pid,
2544 primary,
2545 kTRUE, // check if this is right
2546 kTRUE, // check if this is right
2547 AliAODTrack::kPrimary)
2548 );
2549 aodTrack->ConvertAliPIDtoAODPID();
2550 }
2551 else {
2552 // outside the beam pipe: orphan track
2553 // Don't write them anymore!
2554 continue;
2555 }
2556 } // end of loop on tracks
2557
2558 // muon tracks
2559 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2560 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2561
2562 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2563 p[0] = esdMuTrack->Px();
2564 p[1] = esdMuTrack->Py();
2565 p[2] = esdMuTrack->Pz();
2566 pos[0] = primary->GetX();
2567 pos[1] = primary->GetY();
2568 pos[2] = primary->GetZ();
2569
2570 // has to be changed once the muon pid is provided by the ESD
2571 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2572
2573 primary->AddDaughter(aodTrack =
2574 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2575 0, // no label provided
2576 p,
2577 kTRUE,
2578 pos,
2579 kFALSE,
2580 NULL, // no covariance matrix provided
2581 esdMuTrack->Charge(),
2582 0, // ITSClusterMap is set below
2583 pid,
2584 primary,
2585 kFALSE, // muon tracks are not used to fit the primary vtx
2586 kFALSE, // not used for vertex fit
2587 AliAODTrack::kPrimary)
2588 );
2589
2590 aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
2591 Int_t track2Trigger = esdMuTrack->GetMatchTrigger();
2592 aodTrack->SetMatchTrigger(track2Trigger);
2593 if (track2Trigger)
2594 aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
2595 else
2596 aodTrack->SetChi2MatchTrigger(0.);
2597 }
2598 tracks.Expand(jTracks); // remove 'empty slots' due to unwritten tracks
2599
2600 // Access to the AOD container of PMD clusters
2601 TClonesArray &pmdClusters = *(aod->GetPmdClusters());
2602 Int_t jPmdClusters=0;
2603
2604 for (Int_t iPmd = 0; iPmd < nPmdClus; ++iPmd) {
2605 // file pmd clusters, to be revised!
2606 AliESDPmdTrack *pmdTrack = esd->GetPmdTrack(iPmd);
2607 Int_t nLabel = 0;
2608 Int_t *label = 0x0;
2609 Double_t pos[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ() };
2610 Double_t pid[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; // to be revised!
2611 // type not set!
2612 // assoc cluster not set
2613 new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), pos, pid);
2614 }
2615
2616 // Access to the AOD container of clusters
2617 TClonesArray &caloClusters = *(aod->GetCaloClusters());
2618 Int_t jClusters=0;
2619
2620 // Calo Clusters
2621 TArrayS EMCCellNumber(15000);
2622 TArrayD EMCCellAmplitude(15000);
2623 Int_t nEMCCells = 0;
2624 const Float_t fEMCAmpScale = 1./500;
2625
2626 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
2627
2628 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2629
2630 Int_t id = cluster->GetID();
2631 Int_t nLabel = 0;
2632 Int_t *label = 0x0;
2633 Float_t energy = cluster->E();
2634 cluster->GetPosition(posF);
2635 Char_t ttype=AliAODCluster::kUndef;
2636
2637 if (cluster->GetClusterType() == AliESDCaloCluster::kPHOSCluster) {
2638 ttype=AliAODCluster::kPHOSNeutral;
2639 }
2640 else if (cluster->GetClusterType() == AliESDCaloCluster::kEMCALClusterv1) {
2641 ttype = AliAODCluster::kEMCALClusterv1;
2642 }
2643 else if (cluster->GetClusterType() == AliESDCaloCluster::kEMCALPseudoCluster) {
2644 // Collect raw tower info
2645 for (Int_t iDig = 0; iDig < cluster->GetNumberOfDigits(); iDig++) {
2646 EMCCellNumber[nEMCCells] = cluster->GetDigitIndex()->At(iDig);
2647 EMCCellAmplitude[nEMCCells] = fEMCAmpScale*cluster->GetDigitAmplitude()->At(iDig);
2648 nEMCCells++;
2649 }
2650 // don't write cluster data (it's just a pseudo cluster, holding the tower information)
2651 continue;
2652 }
2653
2654 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
2655 nLabel,
2656 label,
2657 energy,
2658 pos,
2659 NULL,
2660 ttype);
2661
2662 caloCluster->SetCaloCluster(); // to be refined!
2663
2664 }
2665 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
2666 // end of loop on calo clusters
2667
2668 // fill EMC cell info
2669 AliAODCaloCells &EMCCells = *(aod->GetEmcalCells());
2670 EMCCells.CreateContainer(nEMCCells);
2671 EMCCells.SetType(AliAODCaloCells::kEMCAL);
2672 for (Int_t iCell = 0; iCell < nEMCCells; iCell++) {
2673 EMCCells.SetCell(iCell,EMCCellNumber[iCell],EMCCellAmplitude[iCell]);
2674 }
2675 EMCCells.Sort();
2676
2677 // tracklets
2678 AliAODTracklets &SPDTracklets = *(aod->GetTracklets());
2679 const AliMultiplicity *mult = esd->GetMultiplicity();
2680 if (mult) {
2681 if (mult->GetNumberOfTracklets()>0) {
2682 SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
2683
2684 for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
2685 SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n));
2686 }
2687 }
2688 } else {
2689 Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
2690 }
2691
2692 delete [] usedTrack;
2693 delete [] usedV0;
2694 delete [] usedKink;
2695
2696 // fill the tree for this event
2697 aodTree->Fill();
2698 } // end of event loop
2699
2700 aodTree->GetUserInfo()->Add(aod);
2701
2702 // write the tree to the specified file
2703 aodFile = aodTree->GetCurrentFile();
2704 aodFile->cd();
2705 aodTree->Write();
2706
2707 return;
2708}
2709
2710void AliReconstruction::WriteAlignmentData(AliESDEvent* esd)
2711{
2712 // Write space-points which are then used in the alignment procedures
2713 // For the moment only ITS, TRD and TPC
2714
2715 // Load TOF clusters
2716 if (fTracker[3]){
2717 fLoader[3]->LoadRecPoints("read");
2718 TTree* tree = fLoader[3]->TreeR();
2719 if (!tree) {
2720 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2721 return;
2722 }
2723 fTracker[3]->LoadClusters(tree);
2724 }
2725 Int_t ntracks = esd->GetNumberOfTracks();
2726 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2727 {
2728 AliESDtrack *track = esd->GetTrack(itrack);
2729 Int_t nsp = 0;
2730 Int_t idx[200];
2731 for (Int_t iDet = 3; iDet >= 0; iDet--)
2732 nsp += track->GetNcls(iDet);
2733 if (nsp) {
2734 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2735 track->SetTrackPointArray(sp);
2736 Int_t isptrack = 0;
2737 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2738 AliTracker *tracker = fTracker[iDet];
2739 if (!tracker) continue;
2740 Int_t nspdet = track->GetNcls(iDet);
2741 if (nspdet <= 0) continue;
2742 track->GetClusters(iDet,idx);
2743 AliTrackPoint p;
2744 Int_t isp = 0;
2745 Int_t isp2 = 0;
2746 while (isp < nspdet) {
2747 Bool_t isvalid;
2748 if(IsSelected(fgkDetectorName[iDet],fUseTrackingErrorsForAlignment)) {
2749 isvalid = tracker->GetTrackPointTrackingError(idx[isp2],p,track);
2750 } else {
2751 isvalid = tracker->GetTrackPoint(idx[isp2],p);
2752 }
2753 isp2++;
2754 const Int_t kNTPCmax = 159;
2755 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
2756 if (!isvalid) continue;
2757 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2758 }
2759 }
2760 }
2761 }
2762 if (fTracker[3]){
2763 fTracker[3]->UnloadClusters();
2764 fLoader[3]->UnloadRecPoints();
2765 }
2766}
2767
2768//_____________________________________________________________________________
2769void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESDEvent* esd)
2770{
2771 // The method reads the raw-data error log
2772 // accumulated within the rawReader.
2773 // It extracts the raw-data errors related to
2774 // the current event and stores them into
2775 // a TClonesArray inside the esd object.
2776
2777 if (!fRawReader) return;
2778
2779 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2780
2781 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2782 if (!log) continue;
2783 if (iEvent != log->GetEventNumber()) continue;
2784
2785 esd->AddRawDataErrorLog(log);
2786 }
2787
2788}
2789
2790TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2791 // Dump a file content into a char in TNamed
2792 ifstream in;
2793 in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2794 Int_t kBytes = (Int_t)in.tellg();
2795 printf("Size: %d \n",kBytes);
2796 TNamed *fn = 0;
2797 if(in.good()){
2798 char* memblock = new char [kBytes];
2799 in.seekg (0, ios::beg);
2800 in.read (memblock, kBytes);
2801 in.close();
2802 TString fData(memblock,kBytes);
2803 fn = new TNamed(fName,fData);
2804 printf("fData Size: %d \n",fData.Sizeof());
2805 printf("fName Size: %d \n",fName.Sizeof());
2806 printf("fn Size: %d \n",fn->Sizeof());
2807 delete[] memblock;
2808 }
2809 else{
2810 AliInfo(Form("Could not Open %s\n",fPath.Data()));
2811 }
2812
2813 return fn;
2814}
2815
2816void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2817 // This is not really needed in AliReconstruction at the moment
2818 // but can serve as a template
2819
2820 TList *fList = fTree->GetUserInfo();
2821 TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2822 printf("fn Size: %d \n",fn->Sizeof());
2823
2824 TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2825 const char* cdata = fn->GetTitle();
2826 printf("fTmp Size %d\n",fTmp.Sizeof());
2827
2828 int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2829 printf("calculated size %d\n",size);
2830 ofstream out(fName.Data(),ios::out | ios::binary);
2831 out.write(cdata,size);
2832 out.close();
2833
2834}
2835
2836//_____________________________________________________________________________
2837AliQADataMaker * AliReconstruction::GetQADataMaker(Int_t iDet)
2838{
2839 // get the quality assurance data maker object and the loader for a detector
2840
2841 if (fQADataMaker[iDet])
2842 return fQADataMaker[iDet];
2843
2844 // load the QA data maker object
2845 TPluginManager* pluginManager = gROOT->GetPluginManager();
2846 TString detName = fgkDetectorName[iDet];
2847 TString qadmName = "Ali" + detName + "QADataMaker";
2848 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT"))
2849 return NULL;
2850
2851 AliQADataMaker * qadm = NULL;
2852 // first check if a plugin is defined for the quality assurance data maker
2853 TPluginHandler* pluginHandler = pluginManager->FindHandler("AliQADataMaker", detName);
2854 // if not, add a plugin for it
2855 if (!pluginHandler) {
2856 AliDebug(1, Form("defining plugin for %s", qadmName.Data()));
2857 TString libs = gSystem->GetLibraries();
2858 if (libs.Contains("lib" + detName + "base.so") ||
2859 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
2860 pluginManager->AddHandler("AliQADataMaker", detName,
2861 qadmName, detName + "qadm", qadmName + "()");
2862 } else {
2863 pluginManager->AddHandler("AliQADataMaker", detName,
2864 qadmName, detName, qadmName + "()");
2865 }
2866 pluginHandler = pluginManager->FindHandler("AliQADataMaker", detName);
2867 }
2868 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
2869 qadm = (AliQADataMaker *) pluginHandler->ExecPlugin(0);
2870 }
2871 if (qadm) {
2872 AliInfo(Form("Initializing quality assurance data maker for %s", fgkDetectorName[iDet]));
2873 qadm->Init(AliQA::kRECPOINTS, AliCDBManager::Instance()->GetRun(), GetQACycles(fgkDetectorName[iDet]));
2874 qadm->StartOfCycle(AliQA::kRECPOINTS);
2875 qadm->Init(AliQA::kESDS, AliCDBManager::Instance()->GetRun());
2876 qadm->StartOfCycle(AliQA::kESDS, "same") ;
2877 fQADataMaker[iDet] = qadm;
2878 }
2879
2880 return qadm;
2881}
2882
2883//_____________________________________________________________________________
2884Bool_t AliReconstruction::RunQA(const char* detectors, AliESDEvent *& esd)
2885{
2886 // run the Quality Assurance data producer
2887
2888 AliCodeTimerAuto("")
2889 TString detStr = detectors;
2890 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
2891 if (!IsSelected(fgkDetectorName[iDet], detStr))
2892 continue;
2893 AliQADataMaker * qadm = GetQADataMaker(iDet);
2894 if (!qadm)
2895 continue;
2896 AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2897 AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2898
2899 qadm->Exec(AliQA::kESDS, esd) ;
2900 qadm->Increment() ;
2901
2902 AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2903 }
2904 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
2905 AliError(Form("the following detectors were not found: %s",
2906 detStr.Data()));
2907 if (fStopOnError)
2908 return kFALSE;
2909 }
2910
2911 return kTRUE;
2912
2913}
2914
2915//_____________________________________________________________________________
2916void AliReconstruction::CheckQA()
2917{
2918// check the QA of SIM for this run and remove the detectors
2919// with status Fatal
2920
2921 TString newRunLocalReconstruction ;
2922 TString newRunTracking ;
2923 TString newFillESD ;
2924
2925 for (Int_t iDet = 0; iDet < AliQA::kNDET; iDet++) {
2926 TString detName(AliQA::GetDetName(iDet)) ;
2927 AliQA * qa = AliQA::Instance(AliQA::DETECTORINDEX(iDet)) ;
2928 if ( qa->IsSet(AliQA::DETECTORINDEX(iDet), AliQA::kSIM, AliQA::kFATAL)) {
2929 AliInfo(Form("QA status for %s in Hits and/or SDIGITS and/or Digits was Fatal; No reconstruction performed", detName.Data())) ;
2930 } else {
2931 if ( fRunLocalReconstruction.Contains(AliQA::GetDetName(iDet)) ||
2932 fRunLocalReconstruction.Contains("ALL") ) {
2933 newRunLocalReconstruction += detName ;
2934 newRunLocalReconstruction += " " ;
2935 }
2936 if ( fRunTracking.Contains(AliQA::GetDetName(iDet)) ||
2937 fRunTracking.Contains("ALL") ) {
2938 newRunTracking += detName ;
2939 newRunTracking += " " ;
2940 }
2941 if ( fFillESD.Contains(AliQA::GetDetName(iDet)) ||
2942 fFillESD.Contains("ALL") ) {
2943 newFillESD += detName ;
2944 newFillESD += " " ;
2945 }
2946 }
2947 }
2948 fRunLocalReconstruction = newRunLocalReconstruction ;
2949 fRunTracking = newRunTracking ;
2950 fFillESD = newFillESD ;
2951}
2952
2953//_____________________________________________________________________________
2954Int_t AliReconstruction::GetDetIndex(const char* detector)
2955{
2956 // return the detector index corresponding to detector
2957 Int_t index = -1 ;
2958 for (index = 0; index < fgkNDetectors ; index++) {
2959 if ( strcmp(detector, fgkDetectorName[index]) == 0 )
2960 break ;
2961 }
2962 return index ;
2963}