git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Introduce segmented geometry
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
CommitLineData
596a855f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// class for running the reconstruction //
21// //
22// Clusters and tracks are created for all detectors and all events by //
23// typing: //
24// //
25// AliReconstruction rec; //
26// rec.Run(); //
27// //
28// The Run method returns kTRUE in case of successful execution. //
c71de921 29// //
30// If the input to the reconstruction are not simulated digits but raw data, //
31// this can be specified by an argument of the Run method or by the method //
32// //
33// rec.SetInput("..."); //
34// //
35// The input formats and the corresponding argument are: //
36// - DDL raw data files: directory name, ends with "/" //
37// - raw data root file: root file name, extension ".root" //
38// - raw data DATE file: DATE file name, any other non-empty string //
39// - MC root files : empty string, default //
40// //
b26c3770 41// By default all events are reconstructed. The reconstruction can be //
42// limited to a range of events by giving the index of the first and the //
43// last event as an argument to the Run method or by calling //
44// //
45// rec.SetEventRange(..., ...); //
46// //
47// The index -1 (default) can be used for the last event to indicate no //
48// upper limit of the event range. //
49// //
596a855f 50// The name of the galice file can be changed from the default //
e583c30d 51// "galice.root" by passing it as argument to the AliReconstruction //
52// constructor or by //
596a855f 53// //
54// rec.SetGAliceFile("..."); //
55// //
59697224 56// The local reconstruction can be switched on or off for individual //
57// detectors by //
596a855f 58// //
59697224 59// rec.SetRunLocalReconstruction("..."); //
596a855f 60// //
61// The argument is a (case sensitive) string with the names of the //
62// detectors separated by a space. The special string "ALL" selects all //
63// available detectors. This is the default. //
64// //
c71de921 65// The reconstruction of the primary vertex position can be switched off by //
66// //
67// rec.SetRunVertexFinder(kFALSE); //
68// //
b8cd5251 69// The tracking and the creation of ESD tracks can be switched on for //
70// selected detectors by //
596a855f 71// //
b8cd5251 72// rec.SetRunTracking("..."); //
596a855f 73// //
c84a5e9e 74// Uniform/nonuniform field tracking switches (default: uniform field) //
75// //
1d99986f 76// rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
c84a5e9e 77// //
596a855f 78// The filling of additional ESD information can be steered by //
79// //
80// rec.SetFillESD("..."); //
81// //
b8cd5251 82// Again, for both methods the string specifies the list of detectors. //
83// The default is "ALL". //
84// //
85// The call of the shortcut method //
86// //
87// rec.SetRunReconstruction("..."); //
88// //
89// is equivalent to calling SetRunLocalReconstruction, SetRunTracking and //
90// SetFillESD with the same detector selecting string as argument. //
596a855f 91// //
c71de921 92// The reconstruction requires digits or raw data as input. For the creation //
93// of digits and raw data have a look at the class AliSimulation. //
596a855f 94// //
24f7a148 95// For debug purposes the method SetCheckPointLevel can be used. If the //
96// argument is greater than 0, files with ESD events will be written after //
97// selected steps of the reconstruction for each event: //
98// level 1: after tracking and after filling of ESD (final) //
99// level 2: in addition after each tracking step //
100// level 3: in addition after the filling of ESD for each detector //
101// If a final check point file exists for an event, this event will be //
102// skipped in the reconstruction. The tracking and the filling of ESD for //
103// a detector will be skipped as well, if the corresponding check point //
104// file exists. The ESD event will then be loaded from the file instead. //
105// //
596a855f 106///////////////////////////////////////////////////////////////////////////////
107
024a7e64 108#include <TArrayF.h>
109#include <TFile.h>
110#include <TSystem.h>
111#include <TROOT.h>
112#include <TPluginManager.h>
fd46e2d2 113#include <TStopwatch.h>
3103d196 114#include <TGeoManager.h>
2bdb9d38 115#include <TLorentzVector.h>
596a855f 116
117#include "AliReconstruction.h"
b8cd5251 118#include "AliReconstructor.h"
815c2b38 119#include "AliLog.h"
596a855f 120#include "AliRunLoader.h"
121#include "AliRun.h"
b649205a 122#include "AliRawReaderFile.h"
123#include "AliRawReaderDate.h"
124#include "AliRawReaderRoot.h"
596a855f 125#include "AliESD.h"
1d99986f 126#include "AliESDfriend.h"
2257f27e 127#include "AliESDVertex.h"
c84a5e9e 128#include "AliTracker.h"
2257f27e 129#include "AliVertexer.h"
596a855f 130#include "AliHeader.h"
131#include "AliGenEventHeader.h"
b26c3770 132#include "AliPID.h"
596a855f 133#include "AliESDpid.h"
ff8bb5ae 134#include "AliESDtrack.h"
f3a97c86 135
136#include "AliRunTag.h"
f3a97c86 137#include "AliDetectorTag.h"
138#include "AliEventTag.h"
139
98937d93 140#include "AliTrackPointArray.h"
b0314964 141#include "AliCDBManager.h"
6bae477a 142#include "AliCDBEntry.h"
143#include "AliAlignObj.h"
f3a97c86 144
596a855f 145ClassImp(AliReconstruction)
146
147
c757bafd 148//_____________________________________________________________________________
482070f2 149const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
c757bafd 150
596a855f 151//_____________________________________________________________________________
024cf675 152AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 153 const char* name, const char* title) :
154 TNamed(name, title),
155
c84a5e9e 156 fUniformField(kTRUE),
2257f27e 157 fRunVertexFinder(kTRUE),
1f46a9ae 158 fRunHLTTracking(kFALSE),
1d99986f 159 fStopOnError(kFALSE),
160 fWriteAlignmentData(kFALSE),
161 fWriteESDfriend(kFALSE),
162
163 fRunLocalReconstruction("ALL"),
b8cd5251 164 fRunTracking("ALL"),
e583c30d 165 fFillESD("ALL"),
166 fGAliceFileName(gAliceFilename),
b649205a 167 fInput(""),
b26c3770 168 fFirstEvent(0),
169 fLastEvent(-1),
24f7a148 170 fCheckPointLevel(0),
b8cd5251 171 fOptions(),
6bae477a 172 fLoadAlignFromCDB(kTRUE),
173 fLoadAlignData("ALL"),
e583c30d 174
175 fRunLoader(NULL),
b649205a 176 fRawReader(NULL),
b8cd5251 177
98937d93 178 fVertexer(NULL),
179
6bae477a 180 fAlignObjArray(NULL),
024cf675 181 fCDBUri(cdbUri)
596a855f 182{
183// create reconstruction object with default parameters
b8cd5251 184
185 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
186 fReconstructor[iDet] = NULL;
187 fLoader[iDet] = NULL;
188 fTracker[iDet] = NULL;
189 }
e47c4c2e 190 AliPID pid;
596a855f 191}
192
193//_____________________________________________________________________________
194AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 195 TNamed(rec),
196
c84a5e9e 197 fUniformField(rec.fUniformField),
2257f27e 198 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 199 fRunHLTTracking(rec.fRunHLTTracking),
1d99986f 200 fStopOnError(rec.fStopOnError),
201 fWriteAlignmentData(rec.fWriteAlignmentData),
202 fWriteESDfriend(rec.fWriteESDfriend),
203
204 fRunLocalReconstruction(rec.fRunLocalReconstruction),
e583c30d 205 fRunTracking(rec.fRunTracking),
206 fFillESD(rec.fFillESD),
207 fGAliceFileName(rec.fGAliceFileName),
b649205a 208 fInput(rec.fInput),
b26c3770 209 fFirstEvent(rec.fFirstEvent),
210 fLastEvent(rec.fLastEvent),
24f7a148 211 fCheckPointLevel(0),
b8cd5251 212 fOptions(),
6bae477a 213 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
214 fLoadAlignData(rec.fLoadAlignData),
e583c30d 215
216 fRunLoader(NULL),
b649205a 217 fRawReader(NULL),
b8cd5251 218
98937d93 219 fVertexer(NULL),
220
6bae477a 221 fAlignObjArray(rec.fAlignObjArray),
024cf675 222 fCDBUri(rec.fCDBUri)
596a855f 223{
224// copy constructor
225
efd2085e 226 for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
227 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
228 }
b8cd5251 229 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
230 fReconstructor[iDet] = NULL;
231 fLoader[iDet] = NULL;
232 fTracker[iDet] = NULL;
233 }
596a855f 234}
235
236//_____________________________________________________________________________
237AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
238{
239// assignment operator
240
241 this->~AliReconstruction();
242 new(this) AliReconstruction(rec);
243 return *this;
244}
245
246//_____________________________________________________________________________
247AliReconstruction::~AliReconstruction()
248{
249// clean up
250
e583c30d 251 CleanUp();
efd2085e 252 fOptions.Delete();
596a855f 253}
254
024cf675 255//_____________________________________________________________________________
256void AliReconstruction::InitCDBStorage()
257{
258// activate a default CDB storage
259// First check if we have any CDB storage set, because it is used
260// to retrieve the calibration and alignment constants
261
262 AliCDBManager* man = AliCDBManager::Instance();
263 if (!man->IsDefaultStorageSet())
264 {
265 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
266 AliWarning("Default CDB storage not yet set");
267 AliWarning(Form("Using default storage declared in AliSimulation: %s",fCDBUri.Data()));
268 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
269 SetDefaultStorage(fCDBUri);
e08887cd 270
271 Int_t cdbRun = AliCDBManager::Instance()->GetRun();
272 if(cdbRun == -1){
273 AliWarning("AliCDBManager's run number temporarily set to 0!!");
274 AliCDBManager::Instance()->SetRun(0);
275 }
024cf675 276 }
277
278}
279
280//_____________________________________________________________________________
281void AliReconstruction::SetDefaultStorage(const char* uri) {
282// activate a default CDB storage
283
284 AliCDBManager::Instance()->SetDefaultStorage(uri);
285
286}
287
288//_____________________________________________________________________________
289void AliReconstruction::SetSpecificStorage(const char* detName, const char* uri) {
290// activate a detector-specific CDB storage
291
292 AliCDBManager::Instance()->SetSpecificStorage(detName, uri);
293
294}
295
6bae477a 296//_____________________________________________________________________________
297Bool_t AliReconstruction::SetRunNumber()
298{
299 // The method is called in Run() in order
300 // to set a correct run number.
301 // In case of raw data reconstruction the
302 // run number is taken from the raw data header
303
304 if(AliCDBManager::Instance()->GetRun() < 0) {
305 if (!fRunLoader) {
306 AliError("No run loader is found !");
307 return kFALSE;
308 }
309 // read run number from gAlice
310 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
311 AliInfo(Form("Run number: %d",AliCDBManager::Instance()->GetRun()));
312 }
313 return kTRUE;
314}
315
316//_____________________________________________________________________________
317Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
318{
319 // Read collection of alignment objects (AliAlignObj derived) saved
320 // in the TClonesArray ClArrayName and apply them to the geometry
321 // manager singleton.
322 //
323 alObjArray->Sort();
324 Int_t nvols = alObjArray->GetEntriesFast();
325
dc0984f8 326 Bool_t flag = kTRUE;
327
6bae477a 328 for(Int_t j=0; j<nvols; j++)
329 {
330 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
dc0984f8 331 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
6bae477a 332 }
333
334 if (AliDebugLevelClass() >= 1) {
335 gGeoManager->CheckOverlaps(20);
336 TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
337 if(ovexlist->GetEntriesFast()){
338 AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
339 }
340 }
341
dc0984f8 342 return flag;
6bae477a 343
344}
345
346//_____________________________________________________________________________
347Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
348{
349 // Fills array of single detector's alignable objects from CDB
350
351 AliDebug(2, Form("Loading alignment data for detector: %s",detName));
352
353 AliCDBEntry *entry;
354
355 AliCDBPath path(detName,"Align","Data");
356
357 entry=AliCDBManager::Instance()->Get(path.GetPath());
358 if(!entry){
359 AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
360 return kFALSE;
361 }
362 entry->SetOwner(1);
363 TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
364 alignArray->SetOwner(0);
365 AliDebug(2,Form("Found %d alignment objects for %s",
366 alignArray->GetEntries(),detName));
367
368 AliAlignObj *alignObj=0;
369 TIter iter(alignArray);
370
371 // loop over align objects in detector
372 while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
373 fAlignObjArray->Add(alignObj);
374 }
375 // delete entry --- Don't delete, it is cached!
376
377 AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
378 return kTRUE;
379
380}
381
382//_____________________________________________________________________________
383Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
384{
385 // Read the alignment objects from CDB.
386 // Each detector is supposed to have the
387 // alignment objects in DET/Align/Data CDB path.
388 // All the detector objects are then collected,
389 // sorted by geometry level (starting from ALIC) and
390 // then applied to the TGeo geometry.
391 // Finally an overlaps check is performed.
392
393 // Load alignment data from CDB and fill fAlignObjArray
394 if(fLoadAlignFromCDB){
395 if(!fAlignObjArray) fAlignObjArray = new TObjArray();
396
397 //fAlignObjArray->RemoveAll();
398 fAlignObjArray->Clear();
399 fAlignObjArray->SetOwner(0);
400
401 TString detStr = detectors;
402 TString dataNotLoaded="";
403 TString dataLoaded="";
404
405 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
406 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
407 if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
408 dataNotLoaded += fgkDetectorName[iDet];
409 dataNotLoaded += " ";
410 } else {
411 dataLoaded += fgkDetectorName[iDet];
412 dataLoaded += " ";
413 }
414 } // end loop over detectors
415
416 if ((detStr.CompareTo("ALL") == 0)) detStr = "";
417 dataNotLoaded += detStr;
418 AliInfo(Form("Alignment data loaded for: %s",
419 dataLoaded.Data()));
420 AliInfo(Form("Didn't/couldn't load alignment data for: %s",
421 dataNotLoaded.Data()));
422 } // fLoadAlignFromCDB flag
423
424 // Check if the array with alignment objects was
425 // provided by the user. If yes, apply the objects
426 // to the present TGeo geometry
427 if (fAlignObjArray) {
428 if (gGeoManager && gGeoManager->IsClosed()) {
429 if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
dc0984f8 430 AliError("The misalignment of one or more volumes failed!"
431 "Compare the list of simulated detectors and the list of detector alignment data!");
6bae477a 432 return kFALSE;
433 }
434 }
435 else {
436 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
437 return kFALSE;
438 }
439 }
440
441 return kTRUE;
442}
596a855f 443
444//_____________________________________________________________________________
445void AliReconstruction::SetGAliceFile(const char* fileName)
446{
447// set the name of the galice file
448
449 fGAliceFileName = fileName;
450}
451
efd2085e 452//_____________________________________________________________________________
453void AliReconstruction::SetOption(const char* detector, const char* option)
454{
455// set options for the reconstruction of a detector
456
457 TObject* obj = fOptions.FindObject(detector);
458 if (obj) fOptions.Remove(obj);
459 fOptions.Add(new TNamed(detector, option));
460}
461
596a855f 462
463//_____________________________________________________________________________
b26c3770 464Bool_t AliReconstruction::Run(const char* input,
465 Int_t firstEvent, Int_t lastEvent)
596a855f 466{
467// run the reconstruction
468
024cf675 469 InitCDBStorage();
b0314964 470
b649205a 471 // set the input
472 if (!input) input = fInput.Data();
473 TString fileName(input);
474 if (fileName.EndsWith("/")) {
475 fRawReader = new AliRawReaderFile(fileName);
476 } else if (fileName.EndsWith(".root")) {
477 fRawReader = new AliRawReaderRoot(fileName);
478 } else if (!fileName.IsNull()) {
479 fRawReader = new AliRawReaderDate(fileName);
480 fRawReader->SelectEvents(7);
481 }
482
f08fc9f5 483 // get the run loader
484 if (!InitRunLoader()) return kFALSE;
596a855f 485
6bae477a 486 // Set run number in CDBManager (if it is not already set by the user)
487 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
488
489 // Import ideal TGeo geometry and apply misalignment
490 if (!gGeoManager) {
491 TString geom(gSystem->DirName(fGAliceFileName));
492 geom += "/geometry.root";
493 TGeoManager::Import(geom.Data());
494 if (!gGeoManager) if (fStopOnError) return kFALSE;
495 }
496 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
497
498 // Temporary fix by A.Gheata
499 // Could be removed with the next Root version (>5.11)
500 if (gGeoManager) {
501 TIter next(gGeoManager->GetListOfVolumes());
502 TGeoVolume *vol;
503 while ((vol = (TGeoVolume *)next())) {
504 if (vol->GetVoxels()) {
505 if (vol->GetVoxels()->NeedRebuild()) {
506 vol->GetVoxels()->Voxelize();
507 vol->FindOverlaps();
508 }
509 }
510 }
511 }
512
596a855f 513 // local reconstruction
59697224 514 if (!fRunLocalReconstruction.IsNull()) {
515 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 516 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 517 }
518 }
b26c3770 519// if (!fRunVertexFinder && fRunTracking.IsNull() &&
520// fFillESD.IsNull()) return kTRUE;
2257f27e 521
522 // get vertexer
523 if (fRunVertexFinder && !CreateVertexer()) {
524 if (fStopOnError) {
525 CleanUp();
526 return kFALSE;
527 }
528 }
596a855f 529
f08fc9f5 530 // get trackers
b8cd5251 531 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 532 if (fStopOnError) {
533 CleanUp();
534 return kFALSE;
535 }
596a855f 536 }
24f7a148 537
9db3a215 538
539 TStopwatch stopwatch;
540 stopwatch.Start();
541
b26c3770 542 // get the possibly already existing ESD file and tree
1f46a9ae 543 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
b26c3770 544 TFile* fileOld = NULL;
1f46a9ae 545 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 546 if (!gSystem->AccessPathName("AliESDs.root")){
547 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
548 fileOld = TFile::Open("AliESDs.old.root");
549 if (fileOld && fileOld->IsOpen()) {
550 treeOld = (TTree*) fileOld->Get("esdTree");
551 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
1f46a9ae 552 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
553 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
b26c3770 554 }
555 }
556
36711aa4 557 // create the ESD output file and tree
596a855f 558 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
559 if (!file->IsOpen()) {
815c2b38 560 AliError("opening AliESDs.root failed");
b26c3770 561 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 562 }
36711aa4 563 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
564 tree->Branch("ESD", "AliESD", &esd);
1f46a9ae 565 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
566 hlttree->Branch("ESD", "AliESD", &hltesd);
567 delete esd; delete hltesd;
568 esd = NULL; hltesd = NULL;
1d99986f 569
570 // create the file and tree with ESD additions
571 TFile *filef=0; TTree *treef=0; AliESDfriend *esdf=0;
572 if (fWriteESDfriend) {
573 filef = TFile::Open("AliESDfriends.root", "RECREATE");
574 if (!filef->IsOpen()) {
575 AliError("opening AliESDfriends.root failed");
576 }
577 treef = new TTree("esdFriendTree", "Tree with ESD friends");
578 treef->Branch("ESDfriend", "AliESDfriend", &esdf);
579 }
580
36711aa4 581 gROOT->cd();
596a855f 582
583 // loop over events
b649205a 584 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 585
596a855f 586 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 587 if (fRawReader) fRawReader->NextEvent();
588 if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
589 // copy old ESD to the new one
590 if (treeOld) {
591 treeOld->SetBranchAddress("ESD", &esd);
592 treeOld->GetEntry(iEvent);
593 }
594 tree->Fill();
1f46a9ae 595 if (hlttreeOld) {
596 hlttreeOld->SetBranchAddress("ESD", &hltesd);
597 hlttreeOld->GetEntry(iEvent);
598 }
599 hlttree->Fill();
b26c3770 600 continue;
601 }
602
815c2b38 603 AliInfo(Form("processing event %d", iEvent));
596a855f 604 fRunLoader->GetEvent(iEvent);
24f7a148 605
606 char fileName[256];
607 sprintf(fileName, "ESD_%d.%d_final.root",
f08fc9f5 608 fRunLoader->GetHeader()->GetRun(),
609 fRunLoader->GetHeader()->GetEventNrInRun());
24f7a148 610 if (!gSystem->AccessPathName(fileName)) continue;
611
b26c3770 612 // local reconstruction
613 if (!fRunLocalReconstruction.IsNull()) {
614 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
615 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
616 }
617 }
618
1f46a9ae 619 esd = new AliESD; hltesd = new AliESD;
f08fc9f5 620 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 621 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
f08fc9f5 622 esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
1f46a9ae 623 hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
f08fc9f5 624 if (gAlice) {
7ebb06c3 625 esd->SetMagneticField(AliTracker::GetBz());
626 hltesd->SetMagneticField(AliTracker::GetBz());
f08fc9f5 627 } else {
628 // ???
629 }
596a855f 630
2257f27e 631 // vertex finder
632 if (fRunVertexFinder) {
633 if (!ReadESD(esd, "vertex")) {
634 if (!RunVertexFinder(esd)) {
b26c3770 635 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 636 }
637 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
638 }
639 }
640
1f46a9ae 641 // HLT tracking
642 if (!fRunTracking.IsNull()) {
643 if (fRunHLTTracking) {
644 hltesd->SetVertex(esd->GetVertex());
645 if (!RunHLTTracking(hltesd)) {
646 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
647 }
648 }
649 }
650
596a855f 651 // barrel tracking
b8cd5251 652 if (!fRunTracking.IsNull()) {
24f7a148 653 if (!ReadESD(esd, "tracking")) {
654 if (!RunTracking(esd)) {
b26c3770 655 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 656 }
657 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
596a855f 658 }
659 }
660
661 // fill ESD
662 if (!fFillESD.IsNull()) {
663 if (!FillESD(esd, fFillESD)) {
b26c3770 664 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 665 }
666 }
667
668 // combined PID
669 AliESDpid::MakePID(esd);
24f7a148 670 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 671
672 // write ESD
36711aa4 673 tree->Fill();
1f46a9ae 674 // write HLT ESD
675 hlttree->Fill();
24f7a148 676
1d99986f 677 // write ESD friend
678 if (fWriteESDfriend) {
d75007f6 679 esdf=new AliESDfriend();
680 esd->GetESDfriend(esdf);
1d99986f 681 treef->Fill();
682 }
683
f3a97c86 684 if (fCheckPointLevel > 0) WriteESD(esd, "final");
685
1f46a9ae 686 delete esd; delete hltesd;
687 esd = NULL; hltesd = NULL;
596a855f 688 }
689
9db3a215 690 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
691 stopwatch.RealTime(),stopwatch.CpuTime()));
692
36711aa4 693 file->cd();
694 tree->Write();
1f46a9ae 695 hlttree->Write();
f3a97c86 696
1d99986f 697 if (fWriteESDfriend) {
698 filef->cd();
699 treef->Write(); delete treef; filef->Close(); delete filef;
700 }
701
f3a97c86 702 // Create tags for the events in the ESD tree (the ESD tree is always present)
703 // In case of empty events the tags will contain dummy values
704 CreateTag(file);
b26c3770 705 CleanUp(file, fileOld);
596a855f 706
707 return kTRUE;
708}
709
710
711//_____________________________________________________________________________
59697224 712Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 713{
59697224 714// run the local reconstruction
596a855f 715
030b532d 716 TStopwatch stopwatch;
717 stopwatch.Start();
718
596a855f 719 TString detStr = detectors;
b8cd5251 720 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
721 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
722 AliReconstructor* reconstructor = GetReconstructor(iDet);
723 if (!reconstructor) continue;
b26c3770 724 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 725
726 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
727 TStopwatch stopwatchDet;
728 stopwatchDet.Start();
729 if (fRawReader) {
730 fRawReader->RewindEvents();
731 reconstructor->Reconstruct(fRunLoader, fRawReader);
732 } else {
733 reconstructor->Reconstruct(fRunLoader);
596a855f 734 }
5f8272e1 735 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
736 fgkDetectorName[iDet],
737 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
596a855f 738 }
739
740 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 741 AliError(Form("the following detectors were not found: %s",
742 detStr.Data()));
596a855f 743 if (fStopOnError) return kFALSE;
744 }
745
5f8272e1 746 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
747 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 748
596a855f 749 return kTRUE;
750}
751
b26c3770 752//_____________________________________________________________________________
753Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
754{
755// run the local reconstruction
756
757 TStopwatch stopwatch;
758 stopwatch.Start();
759
760 TString detStr = detectors;
761 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
762 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
763 AliReconstructor* reconstructor = GetReconstructor(iDet);
764 if (!reconstructor) continue;
765 AliLoader* loader = fLoader[iDet];
766
767 // conversion of digits
768 if (fRawReader && reconstructor->HasDigitConversion()) {
769 AliInfo(Form("converting raw data digits into root objects for %s",
770 fgkDetectorName[iDet]));
771 TStopwatch stopwatchDet;
772 stopwatchDet.Start();
773 loader->LoadDigits("update");
774 loader->CleanDigits();
775 loader->MakeDigitsContainer();
776 TTree* digitsTree = loader->TreeD();
777 reconstructor->ConvertDigits(fRawReader, digitsTree);
778 loader->WriteDigits("OVERWRITE");
779 loader->UnloadDigits();
5f8272e1 780 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
781 fgkDetectorName[iDet],
782 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 783 }
784
785 // local reconstruction
786 if (!reconstructor->HasLocalReconstruction()) continue;
787 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
788 TStopwatch stopwatchDet;
789 stopwatchDet.Start();
790 loader->LoadRecPoints("update");
791 loader->CleanRecPoints();
792 loader->MakeRecPointsContainer();
793 TTree* clustersTree = loader->TreeR();
794 if (fRawReader && !reconstructor->HasDigitConversion()) {
795 reconstructor->Reconstruct(fRawReader, clustersTree);
796 } else {
797 loader->LoadDigits("read");
798 TTree* digitsTree = loader->TreeD();
799 if (!digitsTree) {
800 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
801 if (fStopOnError) return kFALSE;
802 } else {
803 reconstructor->Reconstruct(digitsTree, clustersTree);
804 }
805 loader->UnloadDigits();
806 }
807 loader->WriteRecPoints("OVERWRITE");
808 loader->UnloadRecPoints();
5f8272e1 809 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
810 fgkDetectorName[iDet],
811 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 812 }
813
814 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815 AliError(Form("the following detectors were not found: %s",
816 detStr.Data()));
817 if (fStopOnError) return kFALSE;
818 }
5f8272e1 819
820 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
821 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 822
823 return kTRUE;
824}
825
596a855f 826//_____________________________________________________________________________
2257f27e 827Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 828{
829// run the barrel tracking
830
030b532d 831 TStopwatch stopwatch;
832 stopwatch.Start();
833
2257f27e 834 AliESDVertex* vertex = NULL;
835 Double_t vtxPos[3] = {0, 0, 0};
836 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
837 TArrayF mcVertex(3);
a6b0b91b 838 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
839 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
840 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
841 }
2257f27e 842
b8cd5251 843 if (fVertexer) {
815c2b38 844 AliInfo("running the ITS vertex finder");
b26c3770 845 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 846 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 847 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 848 if(!vertex){
815c2b38 849 AliWarning("Vertex not found");
c710f220 850 vertex = new AliESDVertex();
d1a50cb5 851 vertex->SetName("default");
2257f27e 852 }
853 else {
854 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 855 vertex->SetName("reconstructed");
2257f27e 856 }
857
858 } else {
815c2b38 859 AliInfo("getting the primary vertex from MC");
2257f27e 860 vertex = new AliESDVertex(vtxPos, vtxErr);
861 }
862
863 if (vertex) {
864 vertex->GetXYZ(vtxPos);
865 vertex->GetSigmaXYZ(vtxErr);
866 } else {
815c2b38 867 AliWarning("no vertex reconstructed");
2257f27e 868 vertex = new AliESDVertex(vtxPos, vtxErr);
869 }
870 esd->SetVertex(vertex);
b8cd5251 871 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
872 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
873 }
2257f27e 874 delete vertex;
875
5f8272e1 876 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
877 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 878
879 return kTRUE;
880}
881
1f46a9ae 882//_____________________________________________________________________________
883Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
884{
885// run the HLT barrel tracking
886
887 TStopwatch stopwatch;
888 stopwatch.Start();
889
890 if (!fRunLoader) {
891 AliError("Missing runLoader!");
892 return kFALSE;
893 }
894
895 AliInfo("running HLT tracking");
896
897 // Get a pointer to the HLT reconstructor
898 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
899 if (!reconstructor) return kFALSE;
900
901 // TPC + ITS
902 for (Int_t iDet = 1; iDet >= 0; iDet--) {
903 TString detName = fgkDetectorName[iDet];
904 AliDebug(1, Form("%s HLT tracking", detName.Data()));
905 reconstructor->SetOption(detName.Data());
906 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
907 if (!tracker) {
908 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
909 if (fStopOnError) return kFALSE;
9dcc06e1 910 continue;
1f46a9ae 911 }
912 Double_t vtxPos[3];
913 Double_t vtxErr[3]={0.005,0.005,0.010};
914 const AliESDVertex *vertex = esd->GetVertex();
915 vertex->GetXYZ(vtxPos);
916 tracker->SetVertex(vtxPos,vtxErr);
917 if(iDet != 1) {
918 fLoader[iDet]->LoadRecPoints("read");
919 TTree* tree = fLoader[iDet]->TreeR();
920 if (!tree) {
921 AliError(Form("Can't get the %s cluster tree", detName.Data()));
922 return kFALSE;
923 }
924 tracker->LoadClusters(tree);
925 }
926 if (tracker->Clusters2Tracks(esd) != 0) {
927 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
928 return kFALSE;
929 }
930 if(iDet != 1) {
931 tracker->UnloadClusters();
932 }
933 delete tracker;
934 }
935
5f8272e1 936 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
937 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 938
939 return kTRUE;
940}
941
2257f27e 942//_____________________________________________________________________________
943Bool_t AliReconstruction::RunTracking(AliESD*& esd)
944{
945// run the barrel tracking
946
947 TStopwatch stopwatch;
948 stopwatch.Start();
24f7a148 949
815c2b38 950 AliInfo("running tracking");
596a855f 951
b8cd5251 952 // pass 1: TPC + ITS inwards
953 for (Int_t iDet = 1; iDet >= 0; iDet--) {
954 if (!fTracker[iDet]) continue;
955 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 956
b8cd5251 957 // load clusters
958 fLoader[iDet]->LoadRecPoints("read");
959 TTree* tree = fLoader[iDet]->TreeR();
960 if (!tree) {
961 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 962 return kFALSE;
963 }
b8cd5251 964 fTracker[iDet]->LoadClusters(tree);
965
966 // run tracking
967 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
968 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 969 return kFALSE;
970 }
b8cd5251 971 if (fCheckPointLevel > 1) {
972 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
973 }
878e1fe1 974 // preliminary PID in TPC needed by the ITS tracker
975 if (iDet == 1) {
976 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 977 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 978 AliESDpid::MakePID(esd);
979 }
b8cd5251 980 }
596a855f 981
b8cd5251 982 // pass 2: ALL backwards
983 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
984 if (!fTracker[iDet]) continue;
985 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
986
987 // load clusters
988 if (iDet > 1) { // all except ITS, TPC
989 TTree* tree = NULL;
7b61cd9c 990 fLoader[iDet]->LoadRecPoints("read");
991 tree = fLoader[iDet]->TreeR();
b8cd5251 992 if (!tree) {
993 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 994 return kFALSE;
995 }
b8cd5251 996 fTracker[iDet]->LoadClusters(tree);
997 }
24f7a148 998
b8cd5251 999 // run tracking
1000 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1001 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1002 return kFALSE;
1003 }
1004 if (fCheckPointLevel > 1) {
1005 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1006 }
24f7a148 1007
b8cd5251 1008 // unload clusters
1009 if (iDet > 2) { // all except ITS, TPC, TRD
1010 fTracker[iDet]->UnloadClusters();
7b61cd9c 1011 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1012 }
8f37df88 1013 // updated PID in TPC needed by the ITS tracker -MI
1014 if (iDet == 1) {
1015 GetReconstructor(1)->FillESD(fRunLoader, esd);
1016 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1017 AliESDpid::MakePID(esd);
1018 }
b8cd5251 1019 }
596a855f 1020
98937d93 1021 // write space-points to the ESD in case alignment data output
1022 // is switched on
1023 if (fWriteAlignmentData)
1024 WriteAlignmentData(esd);
1025
b8cd5251 1026 // pass 3: TRD + TPC + ITS refit inwards
1027 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1028 if (!fTracker[iDet]) continue;
1029 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 1030
b8cd5251 1031 // run tracking
1032 if (fTracker[iDet]->RefitInward(esd) != 0) {
1033 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1034 return kFALSE;
1035 }
1036 if (fCheckPointLevel > 1) {
1037 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1038 }
596a855f 1039
b8cd5251 1040 // unload clusters
1041 fTracker[iDet]->UnloadClusters();
1042 fLoader[iDet]->UnloadRecPoints();
1043 }
ff8bb5ae 1044 //
1045 // Propagate track to the vertex - if not done by ITS
1046 //
1047 Int_t ntracks = esd->GetNumberOfTracks();
1048 for (Int_t itrack=0; itrack<ntracks; itrack++){
1049 const Double_t kRadius = 3; // beam pipe radius
1050 const Double_t kMaxStep = 5; // max step
1051 const Double_t kMaxD = 123456; // max distance to prim vertex
1052 Double_t fieldZ = AliTracker::GetBz(); //
1053 AliESDtrack * track = esd->GetTrack(itrack);
1054 if (!track) continue;
1055 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1056 track->PropagateTo(kRadius, track->GetMass(),kMaxStep,kTRUE);
1057 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1058 }
1059
5f8272e1 1060 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1061 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1062
596a855f 1063 return kTRUE;
1064}
1065
1066//_____________________________________________________________________________
24f7a148 1067Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1068{
1069// fill the event summary data
1070
030b532d 1071 TStopwatch stopwatch;
1072 stopwatch.Start();
815c2b38 1073 AliInfo("filling ESD");
030b532d 1074
596a855f 1075 TString detStr = detectors;
b8cd5251 1076 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1077 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1078 AliReconstructor* reconstructor = GetReconstructor(iDet);
1079 if (!reconstructor) continue;
1080
1081 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1082 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 1083 TTree* clustersTree = NULL;
1084 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1085 fLoader[iDet]->LoadRecPoints("read");
1086 clustersTree = fLoader[iDet]->TreeR();
1087 if (!clustersTree) {
1088 AliError(Form("Can't get the %s clusters tree",
1089 fgkDetectorName[iDet]));
1090 if (fStopOnError) return kFALSE;
1091 }
1092 }
1093 if (fRawReader && !reconstructor->HasDigitConversion()) {
1094 reconstructor->FillESD(fRawReader, clustersTree, esd);
1095 } else {
1096 TTree* digitsTree = NULL;
1097 if (fLoader[iDet]) {
1098 fLoader[iDet]->LoadDigits("read");
1099 digitsTree = fLoader[iDet]->TreeD();
1100 if (!digitsTree) {
1101 AliError(Form("Can't get the %s digits tree",
1102 fgkDetectorName[iDet]));
1103 if (fStopOnError) return kFALSE;
1104 }
1105 }
1106 reconstructor->FillESD(digitsTree, clustersTree, esd);
1107 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1108 }
1109 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1110 fLoader[iDet]->UnloadRecPoints();
1111 }
1112
b8cd5251 1113 if (fRawReader) {
1114 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1115 } else {
1116 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1117 }
b8cd5251 1118 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1119 }
1120 }
1121
1122 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1123 AliError(Form("the following detectors were not found: %s",
1124 detStr.Data()));
596a855f 1125 if (fStopOnError) return kFALSE;
1126 }
1127
5f8272e1 1128 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1129 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1130
596a855f 1131 return kTRUE;
1132}
1133
1134
1135//_____________________________________________________________________________
1136Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1137{
1138// check whether detName is contained in detectors
1139// if yes, it is removed from detectors
1140
1141 // check if all detectors are selected
1142 if ((detectors.CompareTo("ALL") == 0) ||
1143 detectors.BeginsWith("ALL ") ||
1144 detectors.EndsWith(" ALL") ||
1145 detectors.Contains(" ALL ")) {
1146 detectors = "ALL";
1147 return kTRUE;
1148 }
1149
1150 // search for the given detector
1151 Bool_t result = kFALSE;
1152 if ((detectors.CompareTo(detName) == 0) ||
1153 detectors.BeginsWith(detName+" ") ||
1154 detectors.EndsWith(" "+detName) ||
1155 detectors.Contains(" "+detName+" ")) {
1156 detectors.ReplaceAll(detName, "");
1157 result = kTRUE;
1158 }
1159
1160 // clean up the detectors string
1161 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1162 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1163 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1164
1165 return result;
1166}
e583c30d 1167
f08fc9f5 1168//_____________________________________________________________________________
1169Bool_t AliReconstruction::InitRunLoader()
1170{
1171// get or create the run loader
1172
1173 if (gAlice) delete gAlice;
1174 gAlice = NULL;
1175
b26c3770 1176 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1177 // load all base libraries to get the loader classes
1178 TString libs = gSystem->GetLibraries();
1179 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1180 TString detName = fgkDetectorName[iDet];
1181 if (detName == "HLT") continue;
1182 if (libs.Contains("lib" + detName + "base.so")) continue;
1183 gSystem->Load("lib" + detName + "base.so");
1184 }
f08fc9f5 1185 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1186 if (!fRunLoader) {
1187 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1188 CleanUp();
1189 return kFALSE;
1190 }
b26c3770 1191 fRunLoader->CdGAFile();
1192 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1193 if (fRunLoader->LoadgAlice() == 0) {
1194 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1195 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1196 }
f08fc9f5 1197 }
1198 if (!gAlice && !fRawReader) {
1199 AliError(Form("no gAlice object found in file %s",
1200 fGAliceFileName.Data()));
1201 CleanUp();
1202 return kFALSE;
1203 }
1204
1205 } else { // galice.root does not exist
1206 if (!fRawReader) {
1207 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1208 CleanUp();
1209 return kFALSE;
1210 }
1211 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1212 AliConfig::GetDefaultEventFolderName(),
1213 "recreate");
1214 if (!fRunLoader) {
1215 AliError(Form("could not create run loader in file %s",
1216 fGAliceFileName.Data()));
1217 CleanUp();
1218 return kFALSE;
1219 }
1220 fRunLoader->MakeTree("E");
1221 Int_t iEvent = 0;
1222 while (fRawReader->NextEvent()) {
1223 fRunLoader->SetEventNumber(iEvent);
1224 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1225 iEvent, iEvent);
1226 fRunLoader->MakeTree("H");
1227 fRunLoader->TreeE()->Fill();
1228 iEvent++;
1229 }
1230 fRawReader->RewindEvents();
1231 fRunLoader->WriteHeader("OVERWRITE");
1232 fRunLoader->CdGAFile();
1233 fRunLoader->Write(0, TObject::kOverwrite);
1234// AliTracker::SetFieldMap(???);
1235 }
1236
1237 return kTRUE;
1238}
1239
c757bafd 1240//_____________________________________________________________________________
b8cd5251 1241AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1242{
f08fc9f5 1243// get the reconstructor object and the loader for a detector
c757bafd 1244
b8cd5251 1245 if (fReconstructor[iDet]) return fReconstructor[iDet];
1246
1247 // load the reconstructor object
1248 TPluginManager* pluginManager = gROOT->GetPluginManager();
1249 TString detName = fgkDetectorName[iDet];
1250 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1251 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1252
1253 if (detName == "HLT") {
1254 if (!gROOT->GetClass("AliLevel3")) {
1255 gSystem->Load("libAliL3Src.so");
1256 gSystem->Load("libAliL3Misc.so");
1257 gSystem->Load("libAliL3Hough.so");
1258 gSystem->Load("libAliL3Comp.so");
1259 }
1260 }
1261
1262 AliReconstructor* reconstructor = NULL;
1263 // first check if a plugin is defined for the reconstructor
1264 TPluginHandler* pluginHandler =
1265 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1266 // if not, add a plugin for it
1267 if (!pluginHandler) {
b8cd5251 1268 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1269 TString libs = gSystem->GetLibraries();
1270 if (libs.Contains("lib" + detName + "base.so") ||
1271 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1272 pluginManager->AddHandler("AliReconstructor", detName,
1273 recName, detName + "rec", recName + "()");
1274 } else {
1275 pluginManager->AddHandler("AliReconstructor", detName,
1276 recName, detName, recName + "()");
c757bafd 1277 }
b8cd5251 1278 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1279 }
1280 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1281 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1282 }
b8cd5251 1283 if (reconstructor) {
1284 TObject* obj = fOptions.FindObject(detName.Data());
1285 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1286 reconstructor->Init(fRunLoader);
b8cd5251 1287 fReconstructor[iDet] = reconstructor;
1288 }
1289
f08fc9f5 1290 // get or create the loader
1291 if (detName != "HLT") {
1292 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1293 if (!fLoader[iDet]) {
1294 AliConfig::Instance()
1295 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1296 detName, detName);
1297 // first check if a plugin is defined for the loader
1298 TPluginHandler* pluginHandler =
1299 pluginManager->FindHandler("AliLoader", detName);
1300 // if not, add a plugin for it
1301 if (!pluginHandler) {
1302 TString loaderName = "Ali" + detName + "Loader";
1303 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1304 pluginManager->AddHandler("AliLoader", detName,
1305 loaderName, detName + "base",
1306 loaderName + "(const char*, TFolder*)");
1307 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1308 }
1309 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1310 fLoader[iDet] =
1311 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1312 fRunLoader->GetEventFolder());
1313 }
1314 if (!fLoader[iDet]) { // use default loader
1315 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1316 }
1317 if (!fLoader[iDet]) {
1318 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1319 if (fStopOnError) return NULL;
f08fc9f5 1320 } else {
1321 fRunLoader->AddLoader(fLoader[iDet]);
1322 fRunLoader->CdGAFile();
1323 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1324 fRunLoader->Write(0, TObject::kOverwrite);
1325 }
1326 }
1327 }
1328
b8cd5251 1329 return reconstructor;
c757bafd 1330}
1331
2257f27e 1332//_____________________________________________________________________________
1333Bool_t AliReconstruction::CreateVertexer()
1334{
1335// create the vertexer
1336
b8cd5251 1337 fVertexer = NULL;
1338 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1339 if (itsReconstructor) {
b8cd5251 1340 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1341 }
b8cd5251 1342 if (!fVertexer) {
815c2b38 1343 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1344 if (fStopOnError) return kFALSE;
1345 }
1346
1347 return kTRUE;
1348}
1349
24f7a148 1350//_____________________________________________________________________________
b8cd5251 1351Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1352{
f08fc9f5 1353// create the trackers
24f7a148 1354
b8cd5251 1355 TString detStr = detectors;
1356 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1357 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1358 AliReconstructor* reconstructor = GetReconstructor(iDet);
1359 if (!reconstructor) continue;
1360 TString detName = fgkDetectorName[iDet];
1f46a9ae 1361 if (detName == "HLT") {
1362 fRunHLTTracking = kTRUE;
1363 continue;
1364 }
f08fc9f5 1365
1366 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1367 if (!fTracker[iDet] && (iDet < 7)) {
1368 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1369 if (fStopOnError) return kFALSE;
1370 }
1371 }
1372
24f7a148 1373 return kTRUE;
1374}
1375
e583c30d 1376//_____________________________________________________________________________
b26c3770 1377void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1378{
1379// delete trackers and the run loader and close and delete the file
1380
b8cd5251 1381 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1382 delete fReconstructor[iDet];
1383 fReconstructor[iDet] = NULL;
1384 fLoader[iDet] = NULL;
1385 delete fTracker[iDet];
1386 fTracker[iDet] = NULL;
1387 }
1388 delete fVertexer;
1389 fVertexer = NULL;
e583c30d 1390
1391 delete fRunLoader;
1392 fRunLoader = NULL;
b649205a 1393 delete fRawReader;
1394 fRawReader = NULL;
e583c30d 1395
1396 if (file) {
1397 file->Close();
1398 delete file;
1399 }
b26c3770 1400
1401 if (fileOld) {
1402 fileOld->Close();
1403 delete fileOld;
1404 gSystem->Unlink("AliESDs.old.root");
1405 }
e583c30d 1406}
24f7a148 1407
1408
1409//_____________________________________________________________________________
1410Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1411{
1412// read the ESD event from a file
1413
1414 if (!esd) return kFALSE;
1415 char fileName[256];
1416 sprintf(fileName, "ESD_%d.%d_%s.root",
1417 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1418 if (gSystem->AccessPathName(fileName)) return kFALSE;
1419
f3a97c86 1420 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1421 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1422 TFile* file = TFile::Open(fileName);
1423 if (!file || !file->IsOpen()) {
815c2b38 1424 AliError(Form("opening %s failed", fileName));
24f7a148 1425 delete file;
1426 return kFALSE;
1427 }
1428
1429 gROOT->cd();
1430 delete esd;
1431 esd = (AliESD*) file->Get("ESD");
1432 file->Close();
1433 delete file;
1434 return kTRUE;
1435}
1436
1437//_____________________________________________________________________________
1438void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1439{
1440// write the ESD event to a file
1441
1442 if (!esd) return;
1443 char fileName[256];
1444 sprintf(fileName, "ESD_%d.%d_%s.root",
1445 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1446
815c2b38 1447 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1448 TFile* file = TFile::Open(fileName, "recreate");
1449 if (!file || !file->IsOpen()) {
815c2b38 1450 AliError(Form("opening %s failed", fileName));
24f7a148 1451 } else {
1452 esd->Write("ESD");
1453 file->Close();
1454 }
1455 delete file;
1456}
f3a97c86 1457
1458
1459
1460
1461//_____________________________________________________________________________
1462void AliReconstruction::CreateTag(TFile* file)
1463{
2bdb9d38 1464 /////////////
1465 //muon code//
1466 ////////////
56982dd1 1467 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1468 //Variables
56982dd1 1469 Double_t fX,fY,fZ ;
1470 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1471 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1472 Int_t fCharge;
1473 TLorentzVector fEPvector;
1474
1475 Float_t fZVertexCut = 10.0;
1476 Float_t fRhoVertexCut = 2.0;
1477
1478 Float_t fLowPtCut = 1.0;
1479 Float_t fHighPtCut = 3.0;
1480 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1481 ////////////
1482
1483 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1484
089bf903 1485 // Creates the tags for all the events in a given ESD file
4302e20f 1486 Int_t ntrack;
089bf903 1487 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1488 Int_t nPos, nNeg, nNeutr;
1489 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1490 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1491 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1492 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1493 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1494 Int_t fVertexflag;
1495 TString fVertexName;
4302e20f 1496
f3a97c86 1497 AliRunTag *tag = new AliRunTag();
f3a97c86 1498 AliEventTag *evTag = new AliEventTag();
1499 TTree ttag("T","A Tree with event tags");
2bdb9d38 1500 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1501 btag->SetCompressionLevel(9);
2bdb9d38 1502
f3a97c86 1503 AliInfo(Form("Creating the tags......."));
1504
1505 if (!file || !file->IsOpen()) {
1506 AliError(Form("opening failed"));
1507 delete file;
1508 return ;
2bdb9d38 1509 }
1510 Int_t firstEvent = 0,lastEvent = 0;
f3a97c86 1511 TTree *t = (TTree*) file->Get("esdTree");
1512 TBranch * b = t->GetBranch("ESD");
1513 AliESD *esd = 0;
1514 b->SetAddress(&esd);
2bdb9d38 1515
f3a97c86 1516 tag->SetRunId(esd->GetRunNumber());
2bdb9d38 1517
089bf903 1518 Int_t iNumberOfEvents = b->GetEntries();
2bdb9d38 1519 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1520 ntrack = 0;
1521 nPos = 0;
1522 nNeg = 0;
1523 nNeutr =0;
1524 nK0s = 0;
1525 nNeutrons = 0;
1526 nPi0s = 0;
1527 nGamas = 0;
1528 nProtons = 0;
1529 nKaons = 0;
1530 nPions = 0;
1531 nMuons = 0;
1532 nElectrons = 0;
1533 nCh1GeV = 0;
1534 nCh3GeV = 0;
1535 nCh10GeV = 0;
1536 nMu1GeV = 0;
1537 nMu3GeV = 0;
1538 nMu10GeV = 0;
1539 nEl1GeV = 0;
1540 nEl3GeV = 0;
1541 nEl10GeV = 0;
1542 maxPt = .0;
1543 meanPt = .0;
1544 totalP = .0;
d1a50cb5 1545 fVertexflag = 1;
1546
2bdb9d38 1547 b->GetEntry(iEventNumber);
1548 const AliESDVertex * vertexIn = esd->GetVertex();
d1a50cb5 1549 fVertexName = vertexIn->GetName();
1550 if(fVertexName == "default") fVertexflag = 0;
2bdb9d38 1551 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1552 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1553 UInt_t status = esdTrack->GetStatus();
f3a97c86 1554
2bdb9d38 1555 //select only tracks with ITS refit
1556 if ((status&AliESDtrack::kITSrefit)==0) continue;
1557 //select only tracks with TPC refit
1558 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1559
2bdb9d38 1560 //select only tracks with the "combined PID"
1561 if ((status&AliESDtrack::kESDpid)==0) continue;
1562 Double_t p[3];
1563 esdTrack->GetPxPyPz(p);
1564 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1565 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1566 totalP += momentum;
1567 meanPt += fPt;
1568 if(fPt > maxPt) maxPt = fPt;
4302e20f 1569
2bdb9d38 1570 if(esdTrack->GetSign() > 0) {
1571 nPos++;
56982dd1 1572 if(fPt > fLowPtCut) nCh1GeV++;
1573 if(fPt > fHighPtCut) nCh3GeV++;
1574 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1575 }
1576 if(esdTrack->GetSign() < 0) {
1577 nNeg++;
56982dd1 1578 if(fPt > fLowPtCut) nCh1GeV++;
1579 if(fPt > fHighPtCut) nCh3GeV++;
1580 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1581 }
1582 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1583
2bdb9d38 1584 //PID
1585 Double_t prob[5];
1586 esdTrack->GetESDpid(prob);
4302e20f 1587
2bdb9d38 1588 Double_t rcc = 0.0;
1589 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1590 if(rcc == 0.0) continue;
1591 //Bayes' formula
1592 Double_t w[5];
1593 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1594
2bdb9d38 1595 //protons
1596 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1597 //kaons
1598 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1599 //pions
1600 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1601 //electrons
1602 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1603 nElectrons++;
56982dd1 1604 if(fPt > fLowPtCut) nEl1GeV++;
1605 if(fPt > fHighPtCut) nEl3GeV++;
1606 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1607 }
1608 ntrack++;
1609 }//track loop
1610
1611 /////////////
1612 //muon code//
1613 ////////////
1614 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1615 // loop over all reconstructed tracks (also first track of combination)
1616 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1617 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1618 if (muonTrack == 0x0) continue;
4302e20f 1619
2bdb9d38 1620 // Coordinates at vertex
56982dd1 1621 fZ = muonTrack->GetZ();
1622 fY = muonTrack->GetBendingCoor();
1623 fX = muonTrack->GetNonBendingCoor();
4302e20f 1624
56982dd1 1625 fThetaX = muonTrack->GetThetaX();
1626 fThetaY = muonTrack->GetThetaY();
4302e20f 1627
56982dd1 1628 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1629 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1630 fPxRec = fPzRec * TMath::Tan(fThetaX);
1631 fPyRec = fPzRec * TMath::Tan(fThetaY);
1632 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1633
2bdb9d38 1634 //ChiSquare of the track if needed
56982dd1 1635 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1636 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1637 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1638
2bdb9d38 1639 // total number of muons inside a vertex cut
56982dd1 1640 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1641 nMuons++;
56982dd1 1642 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1643 nMu1GeV++;
56982dd1 1644 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1645 nMu3GeV++;
56982dd1 1646 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1647 nMu10GeV++;
1648 }
1649 }
1650 }
1651 }
1652 }//muon track loop
1653
1654 // Fill the event tags
1655 if(ntrack != 0)
1656 meanPt = meanPt/ntrack;
1657
1658 evTag->SetEventId(iEventNumber+1);
1659 evTag->SetVertexX(vertexIn->GetXv());
1660 evTag->SetVertexY(vertexIn->GetYv());
1661 evTag->SetVertexZ(vertexIn->GetZv());
d1a50cb5 1662 evTag->SetVertexZError(vertexIn->GetZRes());
1663 evTag->SetVertexFlag(fVertexflag);
1664
2bdb9d38 1665 evTag->SetT0VertexZ(esd->GetT0zVertex());
1666
1667 evTag->SetTrigger(esd->GetTrigger());
1668
32a5cab4 1669 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1670 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1671 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1672 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1673 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1674 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1675
1676
1677 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1678 evTag->SetNumOfPosTracks(nPos);
1679 evTag->SetNumOfNegTracks(nNeg);
1680 evTag->SetNumOfNeutrTracks(nNeutr);
1681
1682 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1683 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1684 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1685 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1686
1687 evTag->SetNumOfProtons(nProtons);
1688 evTag->SetNumOfKaons(nKaons);
1689 evTag->SetNumOfPions(nPions);
1690 evTag->SetNumOfMuons(nMuons);
1691 evTag->SetNumOfElectrons(nElectrons);
1692 evTag->SetNumOfPhotons(nGamas);
1693 evTag->SetNumOfPi0s(nPi0s);
1694 evTag->SetNumOfNeutrons(nNeutrons);
1695 evTag->SetNumOfKaon0s(nK0s);
1696
1697 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1698 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1699 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1700 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1701 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1702 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1703 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1704 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1705 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1706
85c60a8e 1707 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1708 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1709
1710 evTag->SetTotalMomentum(totalP);
1711 evTag->SetMeanPt(meanPt);
1712 evTag->SetMaxPt(maxPt);
1713
1714 tag->AddEventTag(*evTag);
1715 }
089bf903 1716 lastEvent = iNumberOfEvents;
f3a97c86 1717
1718 ttag.Fill();
1719 tag->Clear();
1720
1721 char fileName[256];
1722 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1723 tag->GetRunId(),firstEvent,lastEvent );
1724 AliInfo(Form("writing tags to file %s", fileName));
1725 AliDebug(1, Form("writing tags to file %s", fileName));
1726
1727 TFile* ftag = TFile::Open(fileName, "recreate");
1728 ftag->cd();
1729 ttag.Write();
1730 ftag->Close();
1731 file->cd();
1732 delete tag;
f3a97c86 1733 delete evTag;
1734}
1735
98937d93 1736void AliReconstruction::WriteAlignmentData(AliESD* esd)
1737{
1738 // Write space-points which are then used in the alignment procedures
1739 // For the moment only ITS, TRD and TPC
1740
1741 // Load TOF clusters
d528ee75 1742 if (fTracker[3]){
1743 fLoader[3]->LoadRecPoints("read");
1744 TTree* tree = fLoader[3]->TreeR();
1745 if (!tree) {
1746 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1747 return;
1748 }
1749 fTracker[3]->LoadClusters(tree);
98937d93 1750 }
98937d93 1751 Int_t ntracks = esd->GetNumberOfTracks();
1752 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1753 {
1754 AliESDtrack *track = esd->GetTrack(itrack);
1755 Int_t nsp = 0;
ef7253ac 1756 Int_t idx[200];
98937d93 1757 for (Int_t iDet = 3; iDet >= 0; iDet--)
1758 nsp += track->GetNcls(iDet);
1759 if (nsp) {
1760 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1761 track->SetTrackPointArray(sp);
1762 Int_t isptrack = 0;
1763 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1764 AliTracker *tracker = fTracker[iDet];
1765 if (!tracker) continue;
1766 Int_t nspdet = track->GetNcls(iDet);
98937d93 1767 if (nspdet <= 0) continue;
1768 track->GetClusters(iDet,idx);
1769 AliTrackPoint p;
1770 Int_t isp = 0;
1771 Int_t isp2 = 0;
1772 while (isp < nspdet) {
1773 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 1774 const Int_t kNTPCmax = 159;
1775 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 1776 if (!isvalid) continue;
1777 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1778 }
98937d93 1779 }
1780 }
1781 }
d528ee75 1782 if (fTracker[3]){
1783 fTracker[3]->UnloadClusters();
1784 fLoader[3]->UnloadRecPoints();
1785 }
98937d93 1786}