]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliReconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
GetEvent() - do not skip the event if the current event is the same, because by defau...
[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) {
679 esdf=new AliESDfriend(*esd);
680 treef->Fill();
681 }
682
f3a97c86 683 if (fCheckPointLevel > 0) WriteESD(esd, "final");
684
1f46a9ae 685 delete esd; delete hltesd;
686 esd = NULL; hltesd = NULL;
596a855f 687 }
688
9db3a215 689 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
690 stopwatch.RealTime(),stopwatch.CpuTime()));
691
36711aa4 692 file->cd();
693 tree->Write();
1f46a9ae 694 hlttree->Write();
f3a97c86 695
1d99986f 696 if (fWriteESDfriend) {
697 filef->cd();
698 treef->Write(); delete treef; filef->Close(); delete filef;
699 }
700
f3a97c86 701 // Create tags for the events in the ESD tree (the ESD tree is always present)
702 // In case of empty events the tags will contain dummy values
703 CreateTag(file);
b26c3770 704 CleanUp(file, fileOld);
596a855f 705
706 return kTRUE;
707}
708
709
710//_____________________________________________________________________________
59697224 711Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 712{
59697224 713// run the local reconstruction
596a855f 714
030b532d 715 TStopwatch stopwatch;
716 stopwatch.Start();
717
596a855f 718 TString detStr = detectors;
b8cd5251 719 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
720 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
721 AliReconstructor* reconstructor = GetReconstructor(iDet);
722 if (!reconstructor) continue;
b26c3770 723 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 724
725 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
726 TStopwatch stopwatchDet;
727 stopwatchDet.Start();
728 if (fRawReader) {
729 fRawReader->RewindEvents();
730 reconstructor->Reconstruct(fRunLoader, fRawReader);
731 } else {
732 reconstructor->Reconstruct(fRunLoader);
596a855f 733 }
5f8272e1 734 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
735 fgkDetectorName[iDet],
736 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
596a855f 737 }
738
739 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 740 AliError(Form("the following detectors were not found: %s",
741 detStr.Data()));
596a855f 742 if (fStopOnError) return kFALSE;
743 }
744
5f8272e1 745 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
746 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 747
596a855f 748 return kTRUE;
749}
750
b26c3770 751//_____________________________________________________________________________
752Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
753{
754// run the local reconstruction
755
756 TStopwatch stopwatch;
757 stopwatch.Start();
758
759 TString detStr = detectors;
760 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
761 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
762 AliReconstructor* reconstructor = GetReconstructor(iDet);
763 if (!reconstructor) continue;
764 AliLoader* loader = fLoader[iDet];
765
766 // conversion of digits
767 if (fRawReader && reconstructor->HasDigitConversion()) {
768 AliInfo(Form("converting raw data digits into root objects for %s",
769 fgkDetectorName[iDet]));
770 TStopwatch stopwatchDet;
771 stopwatchDet.Start();
772 loader->LoadDigits("update");
773 loader->CleanDigits();
774 loader->MakeDigitsContainer();
775 TTree* digitsTree = loader->TreeD();
776 reconstructor->ConvertDigits(fRawReader, digitsTree);
777 loader->WriteDigits("OVERWRITE");
778 loader->UnloadDigits();
5f8272e1 779 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
780 fgkDetectorName[iDet],
781 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 782 }
783
784 // local reconstruction
785 if (!reconstructor->HasLocalReconstruction()) continue;
786 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
787 TStopwatch stopwatchDet;
788 stopwatchDet.Start();
789 loader->LoadRecPoints("update");
790 loader->CleanRecPoints();
791 loader->MakeRecPointsContainer();
792 TTree* clustersTree = loader->TreeR();
793 if (fRawReader && !reconstructor->HasDigitConversion()) {
794 reconstructor->Reconstruct(fRawReader, clustersTree);
795 } else {
796 loader->LoadDigits("read");
797 TTree* digitsTree = loader->TreeD();
798 if (!digitsTree) {
799 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
800 if (fStopOnError) return kFALSE;
801 } else {
802 reconstructor->Reconstruct(digitsTree, clustersTree);
803 }
804 loader->UnloadDigits();
805 }
806 loader->WriteRecPoints("OVERWRITE");
807 loader->UnloadRecPoints();
5f8272e1 808 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
809 fgkDetectorName[iDet],
810 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 811 }
812
813 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
814 AliError(Form("the following detectors were not found: %s",
815 detStr.Data()));
816 if (fStopOnError) return kFALSE;
817 }
5f8272e1 818
819 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
820 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 821
822 return kTRUE;
823}
824
596a855f 825//_____________________________________________________________________________
2257f27e 826Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 827{
828// run the barrel tracking
829
030b532d 830 TStopwatch stopwatch;
831 stopwatch.Start();
832
2257f27e 833 AliESDVertex* vertex = NULL;
834 Double_t vtxPos[3] = {0, 0, 0};
835 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
836 TArrayF mcVertex(3);
a6b0b91b 837 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
838 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
839 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
840 }
2257f27e 841
b8cd5251 842 if (fVertexer) {
815c2b38 843 AliInfo("running the ITS vertex finder");
b26c3770 844 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 845 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 846 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 847 if(!vertex){
815c2b38 848 AliWarning("Vertex not found");
c710f220 849 vertex = new AliESDVertex();
d1a50cb5 850 vertex->SetName("default");
2257f27e 851 }
852 else {
853 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 854 vertex->SetName("reconstructed");
2257f27e 855 }
856
857 } else {
815c2b38 858 AliInfo("getting the primary vertex from MC");
2257f27e 859 vertex = new AliESDVertex(vtxPos, vtxErr);
860 }
861
862 if (vertex) {
863 vertex->GetXYZ(vtxPos);
864 vertex->GetSigmaXYZ(vtxErr);
865 } else {
815c2b38 866 AliWarning("no vertex reconstructed");
2257f27e 867 vertex = new AliESDVertex(vtxPos, vtxErr);
868 }
869 esd->SetVertex(vertex);
b8cd5251 870 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
871 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
872 }
2257f27e 873 delete vertex;
874
5f8272e1 875 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
876 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 877
878 return kTRUE;
879}
880
1f46a9ae 881//_____________________________________________________________________________
882Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
883{
884// run the HLT barrel tracking
885
886 TStopwatch stopwatch;
887 stopwatch.Start();
888
889 if (!fRunLoader) {
890 AliError("Missing runLoader!");
891 return kFALSE;
892 }
893
894 AliInfo("running HLT tracking");
895
896 // Get a pointer to the HLT reconstructor
897 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
898 if (!reconstructor) return kFALSE;
899
900 // TPC + ITS
901 for (Int_t iDet = 1; iDet >= 0; iDet--) {
902 TString detName = fgkDetectorName[iDet];
903 AliDebug(1, Form("%s HLT tracking", detName.Data()));
904 reconstructor->SetOption(detName.Data());
905 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
906 if (!tracker) {
907 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
908 if (fStopOnError) return kFALSE;
9dcc06e1 909 continue;
1f46a9ae 910 }
911 Double_t vtxPos[3];
912 Double_t vtxErr[3]={0.005,0.005,0.010};
913 const AliESDVertex *vertex = esd->GetVertex();
914 vertex->GetXYZ(vtxPos);
915 tracker->SetVertex(vtxPos,vtxErr);
916 if(iDet != 1) {
917 fLoader[iDet]->LoadRecPoints("read");
918 TTree* tree = fLoader[iDet]->TreeR();
919 if (!tree) {
920 AliError(Form("Can't get the %s cluster tree", detName.Data()));
921 return kFALSE;
922 }
923 tracker->LoadClusters(tree);
924 }
925 if (tracker->Clusters2Tracks(esd) != 0) {
926 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
927 return kFALSE;
928 }
929 if(iDet != 1) {
930 tracker->UnloadClusters();
931 }
932 delete tracker;
933 }
934
5f8272e1 935 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
936 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 937
938 return kTRUE;
939}
940
2257f27e 941//_____________________________________________________________________________
942Bool_t AliReconstruction::RunTracking(AliESD*& esd)
943{
944// run the barrel tracking
945
946 TStopwatch stopwatch;
947 stopwatch.Start();
24f7a148 948
815c2b38 949 AliInfo("running tracking");
596a855f 950
b8cd5251 951 // pass 1: TPC + ITS inwards
952 for (Int_t iDet = 1; iDet >= 0; iDet--) {
953 if (!fTracker[iDet]) continue;
954 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 955
b8cd5251 956 // load clusters
957 fLoader[iDet]->LoadRecPoints("read");
958 TTree* tree = fLoader[iDet]->TreeR();
959 if (!tree) {
960 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 961 return kFALSE;
962 }
b8cd5251 963 fTracker[iDet]->LoadClusters(tree);
964
965 // run tracking
966 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
967 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 968 return kFALSE;
969 }
b8cd5251 970 if (fCheckPointLevel > 1) {
971 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
972 }
878e1fe1 973 // preliminary PID in TPC needed by the ITS tracker
974 if (iDet == 1) {
975 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 976 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 977 AliESDpid::MakePID(esd);
978 }
b8cd5251 979 }
596a855f 980
b8cd5251 981 // pass 2: ALL backwards
982 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
983 if (!fTracker[iDet]) continue;
984 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
985
986 // load clusters
987 if (iDet > 1) { // all except ITS, TPC
988 TTree* tree = NULL;
7b61cd9c 989 fLoader[iDet]->LoadRecPoints("read");
990 tree = fLoader[iDet]->TreeR();
b8cd5251 991 if (!tree) {
992 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 993 return kFALSE;
994 }
b8cd5251 995 fTracker[iDet]->LoadClusters(tree);
996 }
24f7a148 997
b8cd5251 998 // run tracking
999 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1000 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1001 return kFALSE;
1002 }
1003 if (fCheckPointLevel > 1) {
1004 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1005 }
24f7a148 1006
b8cd5251 1007 // unload clusters
1008 if (iDet > 2) { // all except ITS, TPC, TRD
1009 fTracker[iDet]->UnloadClusters();
7b61cd9c 1010 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1011 }
8f37df88 1012 // updated PID in TPC needed by the ITS tracker -MI
1013 if (iDet == 1) {
1014 GetReconstructor(1)->FillESD(fRunLoader, esd);
1015 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1016 AliESDpid::MakePID(esd);
1017 }
b8cd5251 1018 }
596a855f 1019
98937d93 1020 // write space-points to the ESD in case alignment data output
1021 // is switched on
1022 if (fWriteAlignmentData)
1023 WriteAlignmentData(esd);
1024
b8cd5251 1025 // pass 3: TRD + TPC + ITS refit inwards
1026 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1027 if (!fTracker[iDet]) continue;
1028 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 1029
b8cd5251 1030 // run tracking
1031 if (fTracker[iDet]->RefitInward(esd) != 0) {
1032 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1033 return kFALSE;
1034 }
1035 if (fCheckPointLevel > 1) {
1036 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1037 }
596a855f 1038
b8cd5251 1039 // unload clusters
1040 fTracker[iDet]->UnloadClusters();
1041 fLoader[iDet]->UnloadRecPoints();
1042 }
ff8bb5ae 1043 //
1044 // Propagate track to the vertex - if not done by ITS
1045 //
1046 Int_t ntracks = esd->GetNumberOfTracks();
1047 for (Int_t itrack=0; itrack<ntracks; itrack++){
1048 const Double_t kRadius = 3; // beam pipe radius
1049 const Double_t kMaxStep = 5; // max step
1050 const Double_t kMaxD = 123456; // max distance to prim vertex
1051 Double_t fieldZ = AliTracker::GetBz(); //
1052 AliESDtrack * track = esd->GetTrack(itrack);
1053 if (!track) continue;
1054 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1055 track->PropagateTo(kRadius, track->GetMass(),kMaxStep,kTRUE);
1056 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1057 }
1058
5f8272e1 1059 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1060 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1061
596a855f 1062 return kTRUE;
1063}
1064
1065//_____________________________________________________________________________
24f7a148 1066Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1067{
1068// fill the event summary data
1069
030b532d 1070 TStopwatch stopwatch;
1071 stopwatch.Start();
815c2b38 1072 AliInfo("filling ESD");
030b532d 1073
596a855f 1074 TString detStr = detectors;
b8cd5251 1075 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1076 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1077 AliReconstructor* reconstructor = GetReconstructor(iDet);
1078 if (!reconstructor) continue;
1079
1080 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1081 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 1082 TTree* clustersTree = NULL;
1083 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1084 fLoader[iDet]->LoadRecPoints("read");
1085 clustersTree = fLoader[iDet]->TreeR();
1086 if (!clustersTree) {
1087 AliError(Form("Can't get the %s clusters tree",
1088 fgkDetectorName[iDet]));
1089 if (fStopOnError) return kFALSE;
1090 }
1091 }
1092 if (fRawReader && !reconstructor->HasDigitConversion()) {
1093 reconstructor->FillESD(fRawReader, clustersTree, esd);
1094 } else {
1095 TTree* digitsTree = NULL;
1096 if (fLoader[iDet]) {
1097 fLoader[iDet]->LoadDigits("read");
1098 digitsTree = fLoader[iDet]->TreeD();
1099 if (!digitsTree) {
1100 AliError(Form("Can't get the %s digits tree",
1101 fgkDetectorName[iDet]));
1102 if (fStopOnError) return kFALSE;
1103 }
1104 }
1105 reconstructor->FillESD(digitsTree, clustersTree, esd);
1106 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1107 }
1108 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1109 fLoader[iDet]->UnloadRecPoints();
1110 }
1111
b8cd5251 1112 if (fRawReader) {
1113 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1114 } else {
1115 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1116 }
b8cd5251 1117 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1118 }
1119 }
1120
1121 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1122 AliError(Form("the following detectors were not found: %s",
1123 detStr.Data()));
596a855f 1124 if (fStopOnError) return kFALSE;
1125 }
1126
5f8272e1 1127 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1128 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1129
596a855f 1130 return kTRUE;
1131}
1132
1133
1134//_____________________________________________________________________________
1135Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1136{
1137// check whether detName is contained in detectors
1138// if yes, it is removed from detectors
1139
1140 // check if all detectors are selected
1141 if ((detectors.CompareTo("ALL") == 0) ||
1142 detectors.BeginsWith("ALL ") ||
1143 detectors.EndsWith(" ALL") ||
1144 detectors.Contains(" ALL ")) {
1145 detectors = "ALL";
1146 return kTRUE;
1147 }
1148
1149 // search for the given detector
1150 Bool_t result = kFALSE;
1151 if ((detectors.CompareTo(detName) == 0) ||
1152 detectors.BeginsWith(detName+" ") ||
1153 detectors.EndsWith(" "+detName) ||
1154 detectors.Contains(" "+detName+" ")) {
1155 detectors.ReplaceAll(detName, "");
1156 result = kTRUE;
1157 }
1158
1159 // clean up the detectors string
1160 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1161 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1162 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1163
1164 return result;
1165}
e583c30d 1166
f08fc9f5 1167//_____________________________________________________________________________
1168Bool_t AliReconstruction::InitRunLoader()
1169{
1170// get or create the run loader
1171
1172 if (gAlice) delete gAlice;
1173 gAlice = NULL;
1174
b26c3770 1175 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1176 // load all base libraries to get the loader classes
1177 TString libs = gSystem->GetLibraries();
1178 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1179 TString detName = fgkDetectorName[iDet];
1180 if (detName == "HLT") continue;
1181 if (libs.Contains("lib" + detName + "base.so")) continue;
1182 gSystem->Load("lib" + detName + "base.so");
1183 }
f08fc9f5 1184 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1185 if (!fRunLoader) {
1186 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1187 CleanUp();
1188 return kFALSE;
1189 }
b26c3770 1190 fRunLoader->CdGAFile();
1191 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1192 if (fRunLoader->LoadgAlice() == 0) {
1193 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1194 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1195 }
f08fc9f5 1196 }
1197 if (!gAlice && !fRawReader) {
1198 AliError(Form("no gAlice object found in file %s",
1199 fGAliceFileName.Data()));
1200 CleanUp();
1201 return kFALSE;
1202 }
1203
1204 } else { // galice.root does not exist
1205 if (!fRawReader) {
1206 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1207 CleanUp();
1208 return kFALSE;
1209 }
1210 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1211 AliConfig::GetDefaultEventFolderName(),
1212 "recreate");
1213 if (!fRunLoader) {
1214 AliError(Form("could not create run loader in file %s",
1215 fGAliceFileName.Data()));
1216 CleanUp();
1217 return kFALSE;
1218 }
1219 fRunLoader->MakeTree("E");
1220 Int_t iEvent = 0;
1221 while (fRawReader->NextEvent()) {
1222 fRunLoader->SetEventNumber(iEvent);
1223 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1224 iEvent, iEvent);
1225 fRunLoader->MakeTree("H");
1226 fRunLoader->TreeE()->Fill();
1227 iEvent++;
1228 }
1229 fRawReader->RewindEvents();
1230 fRunLoader->WriteHeader("OVERWRITE");
1231 fRunLoader->CdGAFile();
1232 fRunLoader->Write(0, TObject::kOverwrite);
1233// AliTracker::SetFieldMap(???);
1234 }
1235
1236 return kTRUE;
1237}
1238
c757bafd 1239//_____________________________________________________________________________
b8cd5251 1240AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1241{
f08fc9f5 1242// get the reconstructor object and the loader for a detector
c757bafd 1243
b8cd5251 1244 if (fReconstructor[iDet]) return fReconstructor[iDet];
1245
1246 // load the reconstructor object
1247 TPluginManager* pluginManager = gROOT->GetPluginManager();
1248 TString detName = fgkDetectorName[iDet];
1249 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1250 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1251
1252 if (detName == "HLT") {
1253 if (!gROOT->GetClass("AliLevel3")) {
1254 gSystem->Load("libAliL3Src.so");
1255 gSystem->Load("libAliL3Misc.so");
1256 gSystem->Load("libAliL3Hough.so");
1257 gSystem->Load("libAliL3Comp.so");
1258 }
1259 }
1260
1261 AliReconstructor* reconstructor = NULL;
1262 // first check if a plugin is defined for the reconstructor
1263 TPluginHandler* pluginHandler =
1264 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1265 // if not, add a plugin for it
1266 if (!pluginHandler) {
b8cd5251 1267 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1268 TString libs = gSystem->GetLibraries();
1269 if (libs.Contains("lib" + detName + "base.so") ||
1270 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1271 pluginManager->AddHandler("AliReconstructor", detName,
1272 recName, detName + "rec", recName + "()");
1273 } else {
1274 pluginManager->AddHandler("AliReconstructor", detName,
1275 recName, detName, recName + "()");
c757bafd 1276 }
b8cd5251 1277 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1278 }
1279 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1280 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1281 }
b8cd5251 1282 if (reconstructor) {
1283 TObject* obj = fOptions.FindObject(detName.Data());
1284 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1285 reconstructor->Init(fRunLoader);
b8cd5251 1286 fReconstructor[iDet] = reconstructor;
1287 }
1288
f08fc9f5 1289 // get or create the loader
1290 if (detName != "HLT") {
1291 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1292 if (!fLoader[iDet]) {
1293 AliConfig::Instance()
1294 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1295 detName, detName);
1296 // first check if a plugin is defined for the loader
1297 TPluginHandler* pluginHandler =
1298 pluginManager->FindHandler("AliLoader", detName);
1299 // if not, add a plugin for it
1300 if (!pluginHandler) {
1301 TString loaderName = "Ali" + detName + "Loader";
1302 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1303 pluginManager->AddHandler("AliLoader", detName,
1304 loaderName, detName + "base",
1305 loaderName + "(const char*, TFolder*)");
1306 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1307 }
1308 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1309 fLoader[iDet] =
1310 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1311 fRunLoader->GetEventFolder());
1312 }
1313 if (!fLoader[iDet]) { // use default loader
1314 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1315 }
1316 if (!fLoader[iDet]) {
1317 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1318 if (fStopOnError) return NULL;
f08fc9f5 1319 } else {
1320 fRunLoader->AddLoader(fLoader[iDet]);
1321 fRunLoader->CdGAFile();
1322 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1323 fRunLoader->Write(0, TObject::kOverwrite);
1324 }
1325 }
1326 }
1327
b8cd5251 1328 return reconstructor;
c757bafd 1329}
1330
2257f27e 1331//_____________________________________________________________________________
1332Bool_t AliReconstruction::CreateVertexer()
1333{
1334// create the vertexer
1335
b8cd5251 1336 fVertexer = NULL;
1337 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1338 if (itsReconstructor) {
b8cd5251 1339 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1340 }
b8cd5251 1341 if (!fVertexer) {
815c2b38 1342 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1343 if (fStopOnError) return kFALSE;
1344 }
1345
1346 return kTRUE;
1347}
1348
24f7a148 1349//_____________________________________________________________________________
b8cd5251 1350Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1351{
f08fc9f5 1352// create the trackers
24f7a148 1353
b8cd5251 1354 TString detStr = detectors;
1355 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1356 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1357 AliReconstructor* reconstructor = GetReconstructor(iDet);
1358 if (!reconstructor) continue;
1359 TString detName = fgkDetectorName[iDet];
1f46a9ae 1360 if (detName == "HLT") {
1361 fRunHLTTracking = kTRUE;
1362 continue;
1363 }
f08fc9f5 1364
1365 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1366 if (!fTracker[iDet] && (iDet < 7)) {
1367 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1368 if (fStopOnError) return kFALSE;
1369 }
1370 }
1371
24f7a148 1372 return kTRUE;
1373}
1374
e583c30d 1375//_____________________________________________________________________________
b26c3770 1376void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1377{
1378// delete trackers and the run loader and close and delete the file
1379
b8cd5251 1380 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1381 delete fReconstructor[iDet];
1382 fReconstructor[iDet] = NULL;
1383 fLoader[iDet] = NULL;
1384 delete fTracker[iDet];
1385 fTracker[iDet] = NULL;
1386 }
1387 delete fVertexer;
1388 fVertexer = NULL;
e583c30d 1389
1390 delete fRunLoader;
1391 fRunLoader = NULL;
b649205a 1392 delete fRawReader;
1393 fRawReader = NULL;
e583c30d 1394
1395 if (file) {
1396 file->Close();
1397 delete file;
1398 }
b26c3770 1399
1400 if (fileOld) {
1401 fileOld->Close();
1402 delete fileOld;
1403 gSystem->Unlink("AliESDs.old.root");
1404 }
e583c30d 1405}
24f7a148 1406
1407
1408//_____________________________________________________________________________
1409Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1410{
1411// read the ESD event from a file
1412
1413 if (!esd) return kFALSE;
1414 char fileName[256];
1415 sprintf(fileName, "ESD_%d.%d_%s.root",
1416 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1417 if (gSystem->AccessPathName(fileName)) return kFALSE;
1418
f3a97c86 1419 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1420 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1421 TFile* file = TFile::Open(fileName);
1422 if (!file || !file->IsOpen()) {
815c2b38 1423 AliError(Form("opening %s failed", fileName));
24f7a148 1424 delete file;
1425 return kFALSE;
1426 }
1427
1428 gROOT->cd();
1429 delete esd;
1430 esd = (AliESD*) file->Get("ESD");
1431 file->Close();
1432 delete file;
1433 return kTRUE;
1434}
1435
1436//_____________________________________________________________________________
1437void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1438{
1439// write the ESD event to a file
1440
1441 if (!esd) return;
1442 char fileName[256];
1443 sprintf(fileName, "ESD_%d.%d_%s.root",
1444 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1445
815c2b38 1446 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1447 TFile* file = TFile::Open(fileName, "recreate");
1448 if (!file || !file->IsOpen()) {
815c2b38 1449 AliError(Form("opening %s failed", fileName));
24f7a148 1450 } else {
1451 esd->Write("ESD");
1452 file->Close();
1453 }
1454 delete file;
1455}
f3a97c86 1456
1457
1458
1459
1460//_____________________________________________________________________________
1461void AliReconstruction::CreateTag(TFile* file)
1462{
2bdb9d38 1463 /////////////
1464 //muon code//
1465 ////////////
56982dd1 1466 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1467 //Variables
56982dd1 1468 Double_t fX,fY,fZ ;
1469 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1470 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1471 Int_t fCharge;
1472 TLorentzVector fEPvector;
1473
1474 Float_t fZVertexCut = 10.0;
1475 Float_t fRhoVertexCut = 2.0;
1476
1477 Float_t fLowPtCut = 1.0;
1478 Float_t fHighPtCut = 3.0;
1479 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1480 ////////////
1481
1482 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1483
089bf903 1484 // Creates the tags for all the events in a given ESD file
4302e20f 1485 Int_t ntrack;
089bf903 1486 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1487 Int_t nPos, nNeg, nNeutr;
1488 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1489 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1490 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1491 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1492 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1493 Int_t fVertexflag;
1494 TString fVertexName;
4302e20f 1495
f3a97c86 1496 AliRunTag *tag = new AliRunTag();
f3a97c86 1497 AliEventTag *evTag = new AliEventTag();
1498 TTree ttag("T","A Tree with event tags");
2bdb9d38 1499 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1500 btag->SetCompressionLevel(9);
2bdb9d38 1501
f3a97c86 1502 AliInfo(Form("Creating the tags......."));
1503
1504 if (!file || !file->IsOpen()) {
1505 AliError(Form("opening failed"));
1506 delete file;
1507 return ;
2bdb9d38 1508 }
1509 Int_t firstEvent = 0,lastEvent = 0;
f3a97c86 1510 TTree *t = (TTree*) file->Get("esdTree");
1511 TBranch * b = t->GetBranch("ESD");
1512 AliESD *esd = 0;
1513 b->SetAddress(&esd);
2bdb9d38 1514
f3a97c86 1515 tag->SetRunId(esd->GetRunNumber());
2bdb9d38 1516
089bf903 1517 Int_t iNumberOfEvents = b->GetEntries();
2bdb9d38 1518 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1519 ntrack = 0;
1520 nPos = 0;
1521 nNeg = 0;
1522 nNeutr =0;
1523 nK0s = 0;
1524 nNeutrons = 0;
1525 nPi0s = 0;
1526 nGamas = 0;
1527 nProtons = 0;
1528 nKaons = 0;
1529 nPions = 0;
1530 nMuons = 0;
1531 nElectrons = 0;
1532 nCh1GeV = 0;
1533 nCh3GeV = 0;
1534 nCh10GeV = 0;
1535 nMu1GeV = 0;
1536 nMu3GeV = 0;
1537 nMu10GeV = 0;
1538 nEl1GeV = 0;
1539 nEl3GeV = 0;
1540 nEl10GeV = 0;
1541 maxPt = .0;
1542 meanPt = .0;
1543 totalP = .0;
d1a50cb5 1544 fVertexflag = 1;
1545
2bdb9d38 1546 b->GetEntry(iEventNumber);
1547 const AliESDVertex * vertexIn = esd->GetVertex();
d1a50cb5 1548 fVertexName = vertexIn->GetName();
1549 if(fVertexName == "default") fVertexflag = 0;
2bdb9d38 1550 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1551 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1552 UInt_t status = esdTrack->GetStatus();
f3a97c86 1553
2bdb9d38 1554 //select only tracks with ITS refit
1555 if ((status&AliESDtrack::kITSrefit)==0) continue;
1556 //select only tracks with TPC refit
1557 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1558
2bdb9d38 1559 //select only tracks with the "combined PID"
1560 if ((status&AliESDtrack::kESDpid)==0) continue;
1561 Double_t p[3];
1562 esdTrack->GetPxPyPz(p);
1563 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1564 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1565 totalP += momentum;
1566 meanPt += fPt;
1567 if(fPt > maxPt) maxPt = fPt;
4302e20f 1568
2bdb9d38 1569 if(esdTrack->GetSign() > 0) {
1570 nPos++;
56982dd1 1571 if(fPt > fLowPtCut) nCh1GeV++;
1572 if(fPt > fHighPtCut) nCh3GeV++;
1573 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1574 }
1575 if(esdTrack->GetSign() < 0) {
1576 nNeg++;
56982dd1 1577 if(fPt > fLowPtCut) nCh1GeV++;
1578 if(fPt > fHighPtCut) nCh3GeV++;
1579 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1580 }
1581 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1582
2bdb9d38 1583 //PID
1584 Double_t prob[5];
1585 esdTrack->GetESDpid(prob);
4302e20f 1586
2bdb9d38 1587 Double_t rcc = 0.0;
1588 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1589 if(rcc == 0.0) continue;
1590 //Bayes' formula
1591 Double_t w[5];
1592 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1593
2bdb9d38 1594 //protons
1595 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1596 //kaons
1597 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1598 //pions
1599 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1600 //electrons
1601 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1602 nElectrons++;
56982dd1 1603 if(fPt > fLowPtCut) nEl1GeV++;
1604 if(fPt > fHighPtCut) nEl3GeV++;
1605 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1606 }
1607 ntrack++;
1608 }//track loop
1609
1610 /////////////
1611 //muon code//
1612 ////////////
1613 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1614 // loop over all reconstructed tracks (also first track of combination)
1615 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1616 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1617 if (muonTrack == 0x0) continue;
4302e20f 1618
2bdb9d38 1619 // Coordinates at vertex
56982dd1 1620 fZ = muonTrack->GetZ();
1621 fY = muonTrack->GetBendingCoor();
1622 fX = muonTrack->GetNonBendingCoor();
4302e20f 1623
56982dd1 1624 fThetaX = muonTrack->GetThetaX();
1625 fThetaY = muonTrack->GetThetaY();
4302e20f 1626
56982dd1 1627 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1628 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1629 fPxRec = fPzRec * TMath::Tan(fThetaX);
1630 fPyRec = fPzRec * TMath::Tan(fThetaY);
1631 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1632
2bdb9d38 1633 //ChiSquare of the track if needed
56982dd1 1634 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1635 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1636 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1637
2bdb9d38 1638 // total number of muons inside a vertex cut
56982dd1 1639 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1640 nMuons++;
56982dd1 1641 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1642 nMu1GeV++;
56982dd1 1643 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1644 nMu3GeV++;
56982dd1 1645 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1646 nMu10GeV++;
1647 }
1648 }
1649 }
1650 }
1651 }//muon track loop
1652
1653 // Fill the event tags
1654 if(ntrack != 0)
1655 meanPt = meanPt/ntrack;
1656
1657 evTag->SetEventId(iEventNumber+1);
1658 evTag->SetVertexX(vertexIn->GetXv());
1659 evTag->SetVertexY(vertexIn->GetYv());
1660 evTag->SetVertexZ(vertexIn->GetZv());
d1a50cb5 1661 evTag->SetVertexZError(vertexIn->GetZRes());
1662 evTag->SetVertexFlag(fVertexflag);
1663
2bdb9d38 1664 evTag->SetT0VertexZ(esd->GetT0zVertex());
1665
1666 evTag->SetTrigger(esd->GetTrigger());
1667
32a5cab4 1668 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1669 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1670 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1671 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1672 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1673 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1674
1675
1676 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1677 evTag->SetNumOfPosTracks(nPos);
1678 evTag->SetNumOfNegTracks(nNeg);
1679 evTag->SetNumOfNeutrTracks(nNeutr);
1680
1681 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1682 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1683 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1684 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1685
1686 evTag->SetNumOfProtons(nProtons);
1687 evTag->SetNumOfKaons(nKaons);
1688 evTag->SetNumOfPions(nPions);
1689 evTag->SetNumOfMuons(nMuons);
1690 evTag->SetNumOfElectrons(nElectrons);
1691 evTag->SetNumOfPhotons(nGamas);
1692 evTag->SetNumOfPi0s(nPi0s);
1693 evTag->SetNumOfNeutrons(nNeutrons);
1694 evTag->SetNumOfKaon0s(nK0s);
1695
1696 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1697 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1698 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1699 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1700 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1701 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1702 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1703 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1704 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1705
85c60a8e 1706 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1707 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1708
1709 evTag->SetTotalMomentum(totalP);
1710 evTag->SetMeanPt(meanPt);
1711 evTag->SetMaxPt(maxPt);
1712
1713 tag->AddEventTag(*evTag);
1714 }
089bf903 1715 lastEvent = iNumberOfEvents;
f3a97c86 1716
1717 ttag.Fill();
1718 tag->Clear();
1719
1720 char fileName[256];
1721 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1722 tag->GetRunId(),firstEvent,lastEvent );
1723 AliInfo(Form("writing tags to file %s", fileName));
1724 AliDebug(1, Form("writing tags to file %s", fileName));
1725
1726 TFile* ftag = TFile::Open(fileName, "recreate");
1727 ftag->cd();
1728 ttag.Write();
1729 ftag->Close();
1730 file->cd();
1731 delete tag;
f3a97c86 1732 delete evTag;
1733}
1734
98937d93 1735void AliReconstruction::WriteAlignmentData(AliESD* esd)
1736{
1737 // Write space-points which are then used in the alignment procedures
1738 // For the moment only ITS, TRD and TPC
1739
1740 // Load TOF clusters
d528ee75 1741 if (fTracker[3]){
1742 fLoader[3]->LoadRecPoints("read");
1743 TTree* tree = fLoader[3]->TreeR();
1744 if (!tree) {
1745 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1746 return;
1747 }
1748 fTracker[3]->LoadClusters(tree);
98937d93 1749 }
98937d93 1750 Int_t ntracks = esd->GetNumberOfTracks();
1751 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1752 {
1753 AliESDtrack *track = esd->GetTrack(itrack);
1754 Int_t nsp = 0;
ef7253ac 1755 Int_t idx[200];
98937d93 1756 for (Int_t iDet = 3; iDet >= 0; iDet--)
1757 nsp += track->GetNcls(iDet);
1758 if (nsp) {
1759 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1760 track->SetTrackPointArray(sp);
1761 Int_t isptrack = 0;
1762 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1763 AliTracker *tracker = fTracker[iDet];
1764 if (!tracker) continue;
1765 Int_t nspdet = track->GetNcls(iDet);
98937d93 1766 if (nspdet <= 0) continue;
1767 track->GetClusters(iDet,idx);
1768 AliTrackPoint p;
1769 Int_t isp = 0;
1770 Int_t isp2 = 0;
1771 while (isp < nspdet) {
1772 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 1773 const Int_t kNTPCmax = 159;
1774 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 1775 if (!isvalid) continue;
1776 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1777 }
98937d93 1778 }
1779 }
1780 }
d528ee75 1781 if (fTracker[3]){
1782 fTracker[3]->UnloadClusters();
1783 fLoader[3]->UnloadRecPoints();
1784 }
98937d93 1785}