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