]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliReconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes to take into account the recent modifications in libESD.pkg
[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// //
973388c2 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// //
596a855f 59// The name of the galice file can be changed from the default //
e583c30d 60// "galice.root" by passing it as argument to the AliReconstruction //
61// constructor or by //
596a855f 62// //
63// rec.SetGAliceFile("..."); //
64// //
59697224 65// The local reconstruction can be switched on or off for individual //
66// detectors by //
596a855f 67// //
59697224 68// rec.SetRunLocalReconstruction("..."); //
596a855f 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// //
c71de921 74// The reconstruction of the primary vertex position can be switched off by //
75// //
76// rec.SetRunVertexFinder(kFALSE); //
77// //
b8cd5251 78// The tracking and the creation of ESD tracks can be switched on for //
79// selected detectors by //
596a855f 80// //
b8cd5251 81// rec.SetRunTracking("..."); //
596a855f 82// //
c84a5e9e 83// Uniform/nonuniform field tracking switches (default: uniform field) //
84// //
1d99986f 85// rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
c84a5e9e 86// //
596a855f 87// The filling of additional ESD information can be steered by //
88// //
89// rec.SetFillESD("..."); //
90// //
b8cd5251 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. //
596a855f 100// //
c71de921 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. //
596a855f 103// //
24f7a148 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// //
596a855f 115///////////////////////////////////////////////////////////////////////////////
116
024a7e64 117#include <TArrayF.h>
118#include <TFile.h>
119#include <TSystem.h>
120#include <TROOT.h>
121#include <TPluginManager.h>
fd46e2d2 122#include <TStopwatch.h>
3103d196 123#include <TGeoManager.h>
2bdb9d38 124#include <TLorentzVector.h>
596a855f 125
126#include "AliReconstruction.h"
b8cd5251 127#include "AliReconstructor.h"
815c2b38 128#include "AliLog.h"
596a855f 129#include "AliRunLoader.h"
130#include "AliRun.h"
b649205a 131#include "AliRawReaderFile.h"
132#include "AliRawReaderDate.h"
133#include "AliRawReaderRoot.h"
001397cd 134#include "AliRawEventHeaderBase.h"
596a855f 135#include "AliESD.h"
1d99986f 136#include "AliESDfriend.h"
2257f27e 137#include "AliESDVertex.h"
32e449be 138#include "AliMultiplicity.h"
c84a5e9e 139#include "AliTracker.h"
2257f27e 140#include "AliVertexer.h"
c5e3e5d1 141#include "AliVertexerTracks.h"
5e4ff34d 142#include "AliV0vertexer.h"
143#include "AliCascadeVertexer.h"
596a855f 144#include "AliHeader.h"
145#include "AliGenEventHeader.h"
b26c3770 146#include "AliPID.h"
596a855f 147#include "AliESDpid.h"
ff8bb5ae 148#include "AliESDtrack.h"
f3a97c86 149
150#include "AliRunTag.h"
f3a97c86 151#include "AliDetectorTag.h"
152#include "AliEventTag.h"
153
98937d93 154#include "AliTrackPointArray.h"
b0314964 155#include "AliCDBManager.h"
6bae477a 156#include "AliCDBEntry.h"
157#include "AliAlignObj.h"
f3a97c86 158
b647652d 159#include "AliCentralTrigger.h"
160#include "AliCTPRawStream.h"
161
596a855f 162ClassImp(AliReconstruction)
163
164
c757bafd 165//_____________________________________________________________________________
b384f8a4 166const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
c757bafd 167
596a855f 168//_____________________________________________________________________________
024cf675 169AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 170 const char* name, const char* title) :
171 TNamed(name, title),
172
c84a5e9e 173 fUniformField(kTRUE),
2257f27e 174 fRunVertexFinder(kTRUE),
1f46a9ae 175 fRunHLTTracking(kFALSE),
1d99986f 176 fStopOnError(kFALSE),
177 fWriteAlignmentData(kFALSE),
178 fWriteESDfriend(kFALSE),
a7807689 179 fWriteAOD(kFALSE),
b647652d 180 fFillTriggerESD(kTRUE),
1d99986f 181
182 fRunLocalReconstruction("ALL"),
b8cd5251 183 fRunTracking("ALL"),
e583c30d 184 fFillESD("ALL"),
185 fGAliceFileName(gAliceFilename),
b649205a 186 fInput(""),
35042093 187 fEquipIdMap(""),
b26c3770 188 fFirstEvent(0),
189 fLastEvent(-1),
973388c2 190 fNumberOfEventsPerFile(1),
24f7a148 191 fCheckPointLevel(0),
b8cd5251 192 fOptions(),
6bae477a 193 fLoadAlignFromCDB(kTRUE),
194 fLoadAlignData("ALL"),
e583c30d 195
196 fRunLoader(NULL),
b649205a 197 fRawReader(NULL),
b8cd5251 198
98937d93 199 fVertexer(NULL),
9178838a 200 fDiamondProfile(NULL),
98937d93 201
6bae477a 202 fAlignObjArray(NULL),
ec92bee0 203 fCDBUri(cdbUri),
204 fSpecCDBUri()
596a855f 205{
206// create reconstruction object with default parameters
b8cd5251 207
208 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
209 fReconstructor[iDet] = NULL;
210 fLoader[iDet] = NULL;
211 fTracker[iDet] = NULL;
212 }
e47c4c2e 213 AliPID pid;
596a855f 214}
215
216//_____________________________________________________________________________
217AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 218 TNamed(rec),
219
c84a5e9e 220 fUniformField(rec.fUniformField),
2257f27e 221 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 222 fRunHLTTracking(rec.fRunHLTTracking),
1d99986f 223 fStopOnError(rec.fStopOnError),
224 fWriteAlignmentData(rec.fWriteAlignmentData),
225 fWriteESDfriend(rec.fWriteESDfriend),
a7807689 226 fWriteAOD(rec.fWriteAOD),
b647652d 227 fFillTriggerESD(rec.fFillTriggerESD),
1d99986f 228
229 fRunLocalReconstruction(rec.fRunLocalReconstruction),
e583c30d 230 fRunTracking(rec.fRunTracking),
231 fFillESD(rec.fFillESD),
232 fGAliceFileName(rec.fGAliceFileName),
b649205a 233 fInput(rec.fInput),
35042093 234 fEquipIdMap(rec.fEquipIdMap),
b26c3770 235 fFirstEvent(rec.fFirstEvent),
236 fLastEvent(rec.fLastEvent),
973388c2 237 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
24f7a148 238 fCheckPointLevel(0),
b8cd5251 239 fOptions(),
6bae477a 240 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
241 fLoadAlignData(rec.fLoadAlignData),
e583c30d 242
243 fRunLoader(NULL),
b649205a 244 fRawReader(NULL),
b8cd5251 245
98937d93 246 fVertexer(NULL),
9178838a 247 fDiamondProfile(NULL),
98937d93 248
6bae477a 249 fAlignObjArray(rec.fAlignObjArray),
ec92bee0 250 fCDBUri(rec.fCDBUri),
251 fSpecCDBUri()
596a855f 252{
253// copy constructor
254
ec92bee0 255 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
efd2085e 256 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
257 }
b8cd5251 258 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
259 fReconstructor[iDet] = NULL;
260 fLoader[iDet] = NULL;
261 fTracker[iDet] = NULL;
262 }
ec92bee0 263 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
264 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
265 }
596a855f 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
e583c30d 283 CleanUp();
efd2085e 284 fOptions.Delete();
ec92bee0 285 fSpecCDBUri.Delete();
596a855f 286}
287
024cf675 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();
ec92bee0 296 if (man->IsDefaultStorageSet())
024cf675 297 {
298 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 299 AliWarning("Default CDB storage has been already set !");
300 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
024cf675 301 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 302 fCDBUri = "";
303 }
304 else {
b8ec52f6 305 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
306 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
307 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 308 man->SetDefaultStorage(fCDBUri);
309 }
310
311 // Now activate the detector specific CDB storage locations
c3a7b59a 312 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
313 TObject* obj = fSpecCDBUri[i];
314 if (!obj) continue;
b8ec52f6 315 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
316 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
317 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
c3a7b59a 318 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
ec92bee0 319 }
727d0766 320 man->Print();
024cf675 321}
322
323//_____________________________________________________________________________
324void AliReconstruction::SetDefaultStorage(const char* uri) {
ec92bee0 325// Store the desired default CDB storage location
326// Activate it later within the Run() method
024cf675 327
ec92bee0 328 fCDBUri = uri;
024cf675 329
330}
331
332//_____________________________________________________________________________
c3a7b59a 333void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
ec92bee0 334// Store a detector-specific CDB storage location
335// Activate it later within the Run() method
024cf675 336
c3a7b59a 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());
ec92bee0 369 if (obj) fSpecCDBUri.Remove(obj);
c3a7b59a 370 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
024cf675 371
372}
373
6bae477a 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
ec92bee0 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()));
6bae477a 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
dc0984f8 421 Bool_t flag = kTRUE;
422
6bae477a 423 for(Int_t j=0; j<nvols; j++)
424 {
425 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
dc0984f8 426 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
6bae477a 427 }
428
429 if (AliDebugLevelClass() >= 1) {
e0625db4 430 gGeoManager->GetTopNode()->CheckOverlaps(20);
6bae477a 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
dc0984f8 437 return flag;
6bae477a 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;
7250c343 513 if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
6bae477a 514 dataLoaded.Data()));
7250c343 515 if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
6bae477a 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) {
dc0984f8 525 AliError("The misalignment of one or more volumes failed!"
526 "Compare the list of simulated detectors and the list of detector alignment data!");
6bae477a 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
8e245d15 536 delete fAlignObjArray; fAlignObjArray=0;
6bae477a 537 return kTRUE;
538}
596a855f 539
540//_____________________________________________________________________________
541void AliReconstruction::SetGAliceFile(const char* fileName)
542{
543// set the name of the galice file
544
545 fGAliceFileName = fileName;
546}
547
efd2085e 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
596a855f 558
559//_____________________________________________________________________________
4a33489c 560Bool_t AliReconstruction::Run(const char* input)
596a855f 561{
562// run the reconstruction
563
b649205a 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 }
35042093 575 if (!fEquipIdMap.IsNull() && fRawReader)
576 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
577
b649205a 578
f08fc9f5 579 // get the run loader
580 if (!InitRunLoader()) return kFALSE;
596a855f 581
ec92bee0 582 // Initialize the CDB storage
583 InitCDBStorage();
584
6bae477a 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 }
8e245d15 595
596 AliCDBManager* man = AliCDBManager::Instance();
6bae477a 597 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
598
596a855f 599 // local reconstruction
59697224 600 if (!fRunLocalReconstruction.IsNull()) {
601 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 602 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 603 }
604 }
b26c3770 605// if (!fRunVertexFinder && fRunTracking.IsNull() &&
606// fFillESD.IsNull()) return kTRUE;
2257f27e 607
608 // get vertexer
609 if (fRunVertexFinder && !CreateVertexer()) {
610 if (fStopOnError) {
611 CleanUp();
612 return kFALSE;
613 }
614 }
596a855f 615
f08fc9f5 616 // get trackers
b8cd5251 617 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 618 if (fStopOnError) {
619 CleanUp();
620 return kFALSE;
621 }
596a855f 622 }
24f7a148 623
9db3a215 624
625 TStopwatch stopwatch;
626 stopwatch.Start();
627
b26c3770 628 // get the possibly already existing ESD file and tree
1f46a9ae 629 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
b26c3770 630 TFile* fileOld = NULL;
1f46a9ae 631 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 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);
1f46a9ae 638 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
639 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
b26c3770 640 }
641 }
642
36711aa4 643 // create the ESD output file and tree
596a855f 644 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
645 if (!file->IsOpen()) {
815c2b38 646 AliError("opening AliESDs.root failed");
b26c3770 647 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 648 }
36711aa4 649 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
650 tree->Branch("ESD", "AliESD", &esd);
1f46a9ae 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;
1d99986f 655
500d54ab 656 // create the branch with ESD additions
657 AliESDfriend *esdf=0;
1d99986f 658 if (fWriteESDfriend) {
500d54ab 659 TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
660 br->SetFile("AliESDfriends.root");
1d99986f 661 }
662
c5e3e5d1 663 AliVertexerTracks tVertexer;
9178838a 664 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
c5e3e5d1 665
596a855f 666 // loop over events
b649205a 667 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 668
596a855f 669 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 670 if (fRawReader) fRawReader->NextEvent();
4a33489c 671 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
b26c3770 672 // copy old ESD to the new one
673 if (treeOld) {
674 treeOld->SetBranchAddress("ESD", &esd);
675 treeOld->GetEntry(iEvent);
676 }
677 tree->Fill();
1f46a9ae 678 if (hlttreeOld) {
679 hlttreeOld->SetBranchAddress("ESD", &hltesd);
680 hlttreeOld->GetEntry(iEvent);
681 }
682 hlttree->Fill();
b26c3770 683 continue;
684 }
685
815c2b38 686 AliInfo(Form("processing event %d", iEvent));
596a855f 687 fRunLoader->GetEvent(iEvent);
24f7a148 688
bb0901a4 689 char aFileName[256];
690 sprintf(aFileName, "ESD_%d.%d_final.root",
f08fc9f5 691 fRunLoader->GetHeader()->GetRun(),
692 fRunLoader->GetHeader()->GetEventNrInRun());
bb0901a4 693 if (!gSystem->AccessPathName(aFileName)) continue;
24f7a148 694
b26c3770 695 // local reconstruction
696 if (!fRunLocalReconstruction.IsNull()) {
697 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
698 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
699 }
700 }
701
1f46a9ae 702 esd = new AliESD; hltesd = new AliESD;
f08fc9f5 703 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 704 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
f08fc9f5 705 esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
1f46a9ae 706 hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
d6ee376f 707
708 // Set magnetic field from the tracker
709 esd->SetMagneticField(AliTracker::GetBz());
710 hltesd->SetMagneticField(AliTracker::GetBz());
596a855f 711
2257f27e 712 // vertex finder
713 if (fRunVertexFinder) {
714 if (!ReadESD(esd, "vertex")) {
715 if (!RunVertexFinder(esd)) {
b26c3770 716 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 717 }
718 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
719 }
720 }
721
1f46a9ae 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
596a855f 732 // barrel tracking
b8cd5251 733 if (!fRunTracking.IsNull()) {
24f7a148 734 if (!ReadESD(esd, "tracking")) {
735 if (!RunTracking(esd)) {
b26c3770 736 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 737 }
738 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
596a855f 739 }
740 }
741
742 // fill ESD
743 if (!fFillESD.IsNull()) {
744 if (!FillESD(esd, fFillESD)) {
b26c3770 745 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 746 }
747 }
001397cd 748 // fill Event header information from the RawEventHeader
749 if (fRawReader){FillRawEventHeaderESD(esd);}
596a855f 750
751 // combined PID
752 AliESDpid::MakePID(esd);
24f7a148 753 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 754
b647652d 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
a03cd2e0 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 }
c5e3e5d1 779
5e4ff34d 780 {
781 // V0 finding
782 AliV0vertexer vtxer;
783 vtxer.Tracks2V0vertices(esd);
784
785 // Cascade finding
786 AliCascadeVertexer cvtxer;
787 cvtxer.V0sTracks2CascadeVertices(esd);
788 }
789
596a855f 790 // write ESD
1d99986f 791 if (fWriteESDfriend) {
d75007f6 792 esdf=new AliESDfriend();
793 esd->GetESDfriend(esdf);
1d99986f 794 }
500d54ab 795 tree->Fill();
796
797 // write HLT ESD
798 hlttree->Fill();
1d99986f 799
f3a97c86 800 if (fCheckPointLevel > 0) WriteESD(esd, "final");
801
500d54ab 802 delete esd; delete esdf; delete hltesd;
803 esd = NULL; esdf=NULL; hltesd = NULL;
596a855f 804 }
805
9db3a215 806 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
807 stopwatch.RealTime(),stopwatch.CpuTime()));
808
36711aa4 809 file->cd();
a9c0e6db 810 if (fWriteESDfriend)
811 tree->SetBranchStatus("ESDfriend*",0);
36711aa4 812 tree->Write();
1f46a9ae 813 hlttree->Write();
f3a97c86 814
a7807689 815 if (fWriteAOD) {
816 CreateAOD(file);
817 }
818
f3a97c86 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);
b26c3770 822 CleanUp(file, fileOld);
596a855f 823
824 return kTRUE;
825}
826
827
828//_____________________________________________________________________________
59697224 829Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 830{
59697224 831// run the local reconstruction
596a855f 832
030b532d 833 TStopwatch stopwatch;
834 stopwatch.Start();
835
8e245d15 836 AliCDBManager* man = AliCDBManager::Instance();
837 Bool_t origCache = man->GetCacheFlag();
838
596a855f 839 TString detStr = detectors;
b8cd5251 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;
b26c3770 844 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 845
846 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
847 TStopwatch stopwatchDet;
848 stopwatchDet.Start();
8e245d15 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
b8cd5251 856 if (fRawReader) {
857 fRawReader->RewindEvents();
858 reconstructor->Reconstruct(fRunLoader, fRawReader);
859 } else {
860 reconstructor->Reconstruct(fRunLoader);
596a855f 861 }
5f8272e1 862 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
863 fgkDetectorName[iDet],
864 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
8e245d15 865
866 // unload calibration data
867 man->ClearCache();
596a855f 868 }
869
8e245d15 870 man->SetCacheFlag(origCache);
871
596a855f 872 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 873 AliError(Form("the following detectors were not found: %s",
874 detStr.Data()));
596a855f 875 if (fStopOnError) return kFALSE;
876 }
877
5f8272e1 878 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
879 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 880
596a855f 881 return kTRUE;
882}
883
b26c3770 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();
5f8272e1 912 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
913 fgkDetectorName[iDet],
914 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 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();
5f8272e1 941 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
942 fgkDetectorName[iDet],
943 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 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 }
5f8272e1 951
952 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
953 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 954
955 return kTRUE;
956}
957
596a855f 958//_____________________________________________________________________________
2257f27e 959Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 960{
961// run the barrel tracking
962
030b532d 963 TStopwatch stopwatch;
964 stopwatch.Start();
965
2257f27e 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);
a6b0b91b 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 }
2257f27e 974
b8cd5251 975 if (fVertexer) {
815c2b38 976 AliInfo("running the ITS vertex finder");
b26c3770 977 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 978 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 979 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 980 if(!vertex){
815c2b38 981 AliWarning("Vertex not found");
c710f220 982 vertex = new AliESDVertex();
d1a50cb5 983 vertex->SetName("default");
2257f27e 984 }
985 else {
986 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 987 vertex->SetName("reconstructed");
2257f27e 988 }
989
990 } else {
815c2b38 991 AliInfo("getting the primary vertex from MC");
2257f27e 992 vertex = new AliESDVertex(vtxPos, vtxErr);
993 }
994
995 if (vertex) {
996 vertex->GetXYZ(vtxPos);
997 vertex->GetSigmaXYZ(vtxErr);
998 } else {
815c2b38 999 AliWarning("no vertex reconstructed");
2257f27e 1000 vertex = new AliESDVertex(vtxPos, vtxErr);
1001 }
1002 esd->SetVertex(vertex);
32e449be 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
b8cd5251 1007 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1008 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1009 }
2257f27e 1010 delete vertex;
1011
5f8272e1 1012 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1013 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 1014
1015 return kTRUE;
1016}
1017
1f46a9ae 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;
9dcc06e1 1046 continue;
1f46a9ae 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
5f8272e1 1072 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1073 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 1074
1075 return kTRUE;
1076}
1077
2257f27e 1078//_____________________________________________________________________________
1079Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1080{
1081// run the barrel tracking
1082
1083 TStopwatch stopwatch;
1084 stopwatch.Start();
24f7a148 1085
815c2b38 1086 AliInfo("running tracking");
596a855f 1087
b8cd5251 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]));
24f7a148 1092
b8cd5251 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]));
24f7a148 1098 return kFALSE;
1099 }
b8cd5251 1100 fTracker[iDet]->LoadClusters(tree);
1101
1102 // run tracking
1103 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1104 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 1105 return kFALSE;
1106 }
b8cd5251 1107 if (fCheckPointLevel > 1) {
1108 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1109 }
878e1fe1 1110 // preliminary PID in TPC needed by the ITS tracker
1111 if (iDet == 1) {
1112 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 1113 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 1114 AliESDpid::MakePID(esd);
1115 }
b8cd5251 1116 }
596a855f 1117
b8cd5251 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;
7b61cd9c 1126 fLoader[iDet]->LoadRecPoints("read");
1127 tree = fLoader[iDet]->TreeR();
b8cd5251 1128 if (!tree) {
1129 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1130 return kFALSE;
1131 }
b8cd5251 1132 fTracker[iDet]->LoadClusters(tree);
1133 }
24f7a148 1134
b8cd5251 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 }
24f7a148 1143
b8cd5251 1144 // unload clusters
1145 if (iDet > 2) { // all except ITS, TPC, TRD
1146 fTracker[iDet]->UnloadClusters();
7b61cd9c 1147 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1148 }
8f37df88 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 }
b8cd5251 1155 }
596a855f 1156
98937d93 1157 // write space-points to the ESD in case alignment data output
1158 // is switched on
1159 if (fWriteAlignmentData)
1160 WriteAlignmentData(esd);
1161
b8cd5251 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]));
596a855f 1166
b8cd5251 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 }
596a855f 1175
b8cd5251 1176 // unload clusters
1177 fTracker[iDet]->UnloadClusters();
1178 fLoader[iDet]->UnloadRecPoints();
1179 }
ff8bb5ae 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;
a28195f7 1192 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
ff8bb5ae 1193 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1194 }
1195
5f8272e1 1196 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1197 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1198
596a855f 1199 return kTRUE;
1200}
1201
1202//_____________________________________________________________________________
24f7a148 1203Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1204{
1205// fill the event summary data
1206
030b532d 1207 TStopwatch stopwatch;
1208 stopwatch.Start();
815c2b38 1209 AliInfo("filling ESD");
030b532d 1210
596a855f 1211 TString detStr = detectors;
b8cd5251 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]));
b26c3770 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
b8cd5251 1249 if (fRawReader) {
1250 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1251 } else {
1252 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1253 }
b8cd5251 1254 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1255 }
1256 }
1257
1258 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1259 AliError(Form("the following detectors were not found: %s",
1260 detStr.Data()));
596a855f 1261 if (fStopOnError) return kFALSE;
1262 }
1263
5f8272e1 1264 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1265 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1266
596a855f 1267 return kTRUE;
1268}
1269
b647652d 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}
596a855f 1309
001397cd 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
596a855f 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}
e583c30d 1366
f08fc9f5 1367//_____________________________________________________________________________
1368Bool_t AliReconstruction::InitRunLoader()
1369{
1370// get or create the run loader
1371
1372 if (gAlice) delete gAlice;
1373 gAlice = NULL;
1374
b26c3770 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 }
f08fc9f5 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 }
b26c3770 1390 fRunLoader->CdGAFile();
1391 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1392 if (fRunLoader->LoadgAlice() == 0) {
1393 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1394 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1395 }
f08fc9f5 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();
973388c2 1430 if (fNumberOfEventsPerFile > 0)
1431 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1432 else
1433 fRunLoader->SetNumberOfEventsPerFile(iEvent);
f08fc9f5 1434 fRunLoader->WriteHeader("OVERWRITE");
1435 fRunLoader->CdGAFile();
1436 fRunLoader->Write(0, TObject::kOverwrite);
1437// AliTracker::SetFieldMap(???);
1438 }
1439
1440 return kTRUE;
1441}
1442
c757bafd 1443//_____________________________________________________________________________
b8cd5251 1444AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1445{
f08fc9f5 1446// get the reconstructor object and the loader for a detector
c757bafd 1447
b8cd5251 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";
f08fc9f5 1454 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1455
1456 if (detName == "HLT") {
1457 if (!gROOT->GetClass("AliLevel3")) {
4aa41877 1458 gSystem->Load("libAliHLTSrc.so");
1459 gSystem->Load("libAliHLTMisc.so");
1460 gSystem->Load("libAliHLTHough.so");
1461 gSystem->Load("libAliHLTComp.so");
b8cd5251 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);
f08fc9f5 1469 // if not, add a plugin for it
1470 if (!pluginHandler) {
b8cd5251 1471 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1472 TString libs = gSystem->GetLibraries();
1473 if (libs.Contains("lib" + detName + "base.so") ||
1474 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1475 pluginManager->AddHandler("AliReconstructor", detName,
1476 recName, detName + "rec", recName + "()");
1477 } else {
1478 pluginManager->AddHandler("AliReconstructor", detName,
1479 recName, detName, recName + "()");
c757bafd 1480 }
b8cd5251 1481 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1482 }
1483 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1484 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1485 }
b8cd5251 1486 if (reconstructor) {
1487 TObject* obj = fOptions.FindObject(detName.Data());
1488 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1489 reconstructor->Init(fRunLoader);
b8cd5251 1490 fReconstructor[iDet] = reconstructor;
1491 }
1492
f08fc9f5 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
bb0901a4 1501 pluginHandler =
f08fc9f5 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()));
6667b602 1522 if (fStopOnError) return NULL;
f08fc9f5 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
b8cd5251 1532 return reconstructor;
c757bafd 1533}
1534
2257f27e 1535//_____________________________________________________________________________
1536Bool_t AliReconstruction::CreateVertexer()
1537{
1538// create the vertexer
1539
b8cd5251 1540 fVertexer = NULL;
1541 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1542 if (itsReconstructor) {
b8cd5251 1543 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1544 }
b8cd5251 1545 if (!fVertexer) {
815c2b38 1546 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1547 if (fStopOnError) return kFALSE;
1548 }
1549
1550 return kTRUE;
1551}
1552
24f7a148 1553//_____________________________________________________________________________
b8cd5251 1554Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1555{
f08fc9f5 1556// create the trackers
24f7a148 1557
b8cd5251 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];
1f46a9ae 1564 if (detName == "HLT") {
1565 fRunHLTTracking = kTRUE;
1566 continue;
1567 }
f08fc9f5 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()));
8250d5f5 1572 if (fStopOnError) return kFALSE;
1573 }
1574 }
1575
24f7a148 1576 return kTRUE;
1577}
1578
e583c30d 1579//_____________________________________________________________________________
b26c3770 1580void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1581{
1582// delete trackers and the run loader and close and delete the file
1583
b8cd5251 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;
9178838a 1593 delete fDiamondProfile;
1594 fDiamondProfile = NULL;
e583c30d 1595
1596 delete fRunLoader;
1597 fRunLoader = NULL;
b649205a 1598 delete fRawReader;
1599 fRawReader = NULL;
e583c30d 1600
1601 if (file) {
1602 file->Close();
1603 delete file;
1604 }
b26c3770 1605
1606 if (fileOld) {
1607 fileOld->Close();
1608 delete fileOld;
1609 gSystem->Unlink("AliESDs.old.root");
1610 }
e583c30d 1611}
24f7a148 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
f3a97c86 1625 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1626 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1627 TFile* file = TFile::Open(fileName);
1628 if (!file || !file->IsOpen()) {
815c2b38 1629 AliError(Form("opening %s failed", fileName));
24f7a148 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
815c2b38 1652 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1653 TFile* file = TFile::Open(fileName, "recreate");
1654 if (!file || !file->IsOpen()) {
815c2b38 1655 AliError(Form("opening %s failed", fileName));
24f7a148 1656 } else {
1657 esd->Write("ESD");
1658 file->Close();
1659 }
1660 delete file;
1661}
f3a97c86 1662
1663
1664
1665
1666//_____________________________________________________________________________
1667void AliReconstruction::CreateTag(TFile* file)
1668{
2bdb9d38 1669 /////////////
1670 //muon code//
1671 ////////////
56982dd1 1672 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1673 //Variables
56982dd1 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;
2bdb9d38 1686 ////////////
1687
1688 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1689
089bf903 1690 // Creates the tags for all the events in a given ESD file
4302e20f 1691 Int_t ntrack;
089bf903 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;
d1a50cb5 1699 Int_t fVertexflag;
1387d165 1700 Int_t iRunNumber = 0;
547d75a6 1701 TString fVertexName("default");
4302e20f 1702
f3a97c86 1703 AliRunTag *tag = new AliRunTag();
f3a97c86 1704 AliEventTag *evTag = new AliEventTag();
1705 TTree ttag("T","A Tree with event tags");
2bdb9d38 1706 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1707 btag->SetCompressionLevel(9);
2bdb9d38 1708
f3a97c86 1709 AliInfo(Form("Creating the tags......."));
1710
1711 if (!file || !file->IsOpen()) {
1712 AliError(Form("opening failed"));
1713 delete file;
1714 return ;
2bdb9d38 1715 }
d71aa170 1716 Int_t lastEvent = 0;
f3a97c86 1717 TTree *t = (TTree*) file->Get("esdTree");
1718 TBranch * b = t->GetBranch("ESD");
1719 AliESD *esd = 0;
1720 b->SetAddress(&esd);
1387d165 1721
d71aa170 1722 b->GetEntry(fFirstEvent);
1387d165 1723 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1724
bb0901a4 1725 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1726 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1727 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 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;
547d75a6 1753 fVertexflag = 0;
d1a50cb5 1754
2bdb9d38 1755 b->GetEntry(iEventNumber);
1387d165 1756 iRunNumber = esd->GetRunNumber();
1757 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1758 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1759 if (!vertexIn) AliError("ESD has not defined vertex.");
1760 if (vertexIn) fVertexName = vertexIn->GetName();
1761 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1762 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1763 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1764 UInt_t status = esdTrack->GetStatus();
f3a97c86 1765
2bdb9d38 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;
4302e20f 1770
2bdb9d38 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;
4302e20f 1780
2bdb9d38 1781 if(esdTrack->GetSign() > 0) {
1782 nPos++;
56982dd1 1783 if(fPt > fLowPtCut) nCh1GeV++;
1784 if(fPt > fHighPtCut) nCh3GeV++;
1785 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1786 }
1787 if(esdTrack->GetSign() < 0) {
1788 nNeg++;
56982dd1 1789 if(fPt > fLowPtCut) nCh1GeV++;
1790 if(fPt > fHighPtCut) nCh3GeV++;
1791 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1792 }
1793 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1794
2bdb9d38 1795 //PID
1796 Double_t prob[5];
1797 esdTrack->GetESDpid(prob);
4302e20f 1798
2bdb9d38 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;
4302e20f 1805
2bdb9d38 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++;
56982dd1 1815 if(fPt > fLowPtCut) nEl1GeV++;
1816 if(fPt > fHighPtCut) nEl3GeV++;
1817 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 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;
4302e20f 1830
2bdb9d38 1831 // Coordinates at vertex
56982dd1 1832 fZ = muonTrack->GetZ();
1833 fY = muonTrack->GetBendingCoor();
1834 fX = muonTrack->GetNonBendingCoor();
4302e20f 1835
56982dd1 1836 fThetaX = muonTrack->GetThetaX();
1837 fThetaY = muonTrack->GetThetaY();
4302e20f 1838
56982dd1 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()));
f3a97c86 1844
2bdb9d38 1845 //ChiSquare of the track if needed
56982dd1 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);
4302e20f 1849
2bdb9d38 1850 // total number of muons inside a vertex cut
56982dd1 1851 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1852 nMuons++;
56982dd1 1853 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1854 nMu1GeV++;
56982dd1 1855 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1856 nMu3GeV++;
56982dd1 1857 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 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);
547d75a6 1870 if (vertexIn) {
1871 evTag->SetVertexX(vertexIn->GetXv());
1872 evTag->SetVertexY(vertexIn->GetYv());
1873 evTag->SetVertexZ(vertexIn->GetZv());
1874 evTag->SetVertexZError(vertexIn->GetZRes());
1875 }
d1a50cb5 1876 evTag->SetVertexFlag(fVertexflag);
1877
2bdb9d38 1878 evTag->SetT0VertexZ(esd->GetT0zVertex());
1879
8bd8ac26 1880 evTag->SetTriggerMask(esd->GetTriggerMask());
1881 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1882
32a5cab4 1883 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1884 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1885 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1886 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 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
85c60a8e 1921 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1922 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1923
1924 evTag->SetTotalMomentum(totalP);
1925 evTag->SetMeanPt(meanPt);
1926 evTag->SetMaxPt(maxPt);
1927
1387d165 1928 tag->SetRunId(iInitRunNumber);
2bdb9d38 1929 tag->AddEventTag(*evTag);
1930 }
bb0901a4 1931 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 1932 else lastEvent = fLastEvent;
f3a97c86 1933
1934 ttag.Fill();
1935 tag->Clear();
1936
1937 char fileName[256];
1938 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 1939 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 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;
f3a97c86 1949 delete evTag;
1950}
1951
a7807689 1952//_____________________________________________________________________________
1953void AliReconstruction::CreateAOD(TFile* esdFile)
1954{
1955 // do nothing for now
1956
1957 return;
1958}
1959
1960
98937d93 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
d528ee75 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);
98937d93 1975 }
98937d93 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;
ef7253ac 1981 Int_t idx[200];
98937d93 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);
98937d93 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++;
160db090 1999 const Int_t kNTPCmax = 159;
2000 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2001 if (!isvalid) continue;
2002 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2003 }
98937d93 2004 }
2005 }
2006 }
d528ee75 2007 if (fTracker[3]){
2008 fTracker[3]->UnloadClusters();
2009 fLoader[3]->UnloadRecPoints();
2010 }
98937d93 2011}