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