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