]> 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
Compilation of AliAnalysisGoodies.cxx only if Root was compiled with XML support
[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());
31fd97b2 705 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
706 hltesd->SetEventNumberInFile(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");
31fd97b2 1322 esd->SetBunchCrossNumber(0);
1323 esd->SetOrbitNumber(0);
001397cd 1324 esd->SetTimeStamp(0);
1325 esd->SetEventType(0);
1326 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1327 if (eventHeader){
31fd97b2 1328 esd->SetBunchCrossNumber((eventHeader->GetP("Id")[0]));
1329 esd->SetOrbitNumber((eventHeader->GetP("Id")[1]));
001397cd 1330 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
31fd97b2 1331 esd->SetEventType((eventHeader->Get("Type")));
001397cd 1332 }
1333
1334 return kTRUE;
1335}
1336
1337
596a855f 1338//_____________________________________________________________________________
1339Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1340{
1341// check whether detName is contained in detectors
1342// if yes, it is removed from detectors
1343
1344 // check if all detectors are selected
1345 if ((detectors.CompareTo("ALL") == 0) ||
1346 detectors.BeginsWith("ALL ") ||
1347 detectors.EndsWith(" ALL") ||
1348 detectors.Contains(" ALL ")) {
1349 detectors = "ALL";
1350 return kTRUE;
1351 }
1352
1353 // search for the given detector
1354 Bool_t result = kFALSE;
1355 if ((detectors.CompareTo(detName) == 0) ||
1356 detectors.BeginsWith(detName+" ") ||
1357 detectors.EndsWith(" "+detName) ||
1358 detectors.Contains(" "+detName+" ")) {
1359 detectors.ReplaceAll(detName, "");
1360 result = kTRUE;
1361 }
1362
1363 // clean up the detectors string
1364 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1365 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1366 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1367
1368 return result;
1369}
e583c30d 1370
f08fc9f5 1371//_____________________________________________________________________________
1372Bool_t AliReconstruction::InitRunLoader()
1373{
1374// get or create the run loader
1375
1376 if (gAlice) delete gAlice;
1377 gAlice = NULL;
1378
b26c3770 1379 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1380 // load all base libraries to get the loader classes
1381 TString libs = gSystem->GetLibraries();
1382 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1383 TString detName = fgkDetectorName[iDet];
1384 if (detName == "HLT") continue;
1385 if (libs.Contains("lib" + detName + "base.so")) continue;
1386 gSystem->Load("lib" + detName + "base.so");
1387 }
f08fc9f5 1388 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1389 if (!fRunLoader) {
1390 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1391 CleanUp();
1392 return kFALSE;
1393 }
b26c3770 1394 fRunLoader->CdGAFile();
1395 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1396 if (fRunLoader->LoadgAlice() == 0) {
1397 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1398 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1399 }
f08fc9f5 1400 }
1401 if (!gAlice && !fRawReader) {
1402 AliError(Form("no gAlice object found in file %s",
1403 fGAliceFileName.Data()));
1404 CleanUp();
1405 return kFALSE;
1406 }
1407
1408 } else { // galice.root does not exist
1409 if (!fRawReader) {
1410 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1411 CleanUp();
1412 return kFALSE;
1413 }
1414 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1415 AliConfig::GetDefaultEventFolderName(),
1416 "recreate");
1417 if (!fRunLoader) {
1418 AliError(Form("could not create run loader in file %s",
1419 fGAliceFileName.Data()));
1420 CleanUp();
1421 return kFALSE;
1422 }
1423 fRunLoader->MakeTree("E");
1424 Int_t iEvent = 0;
1425 while (fRawReader->NextEvent()) {
1426 fRunLoader->SetEventNumber(iEvent);
1427 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1428 iEvent, iEvent);
1429 fRunLoader->MakeTree("H");
1430 fRunLoader->TreeE()->Fill();
1431 iEvent++;
1432 }
1433 fRawReader->RewindEvents();
973388c2 1434 if (fNumberOfEventsPerFile > 0)
1435 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1436 else
1437 fRunLoader->SetNumberOfEventsPerFile(iEvent);
f08fc9f5 1438 fRunLoader->WriteHeader("OVERWRITE");
1439 fRunLoader->CdGAFile();
1440 fRunLoader->Write(0, TObject::kOverwrite);
1441// AliTracker::SetFieldMap(???);
1442 }
1443
1444 return kTRUE;
1445}
1446
c757bafd 1447//_____________________________________________________________________________
b8cd5251 1448AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1449{
f08fc9f5 1450// get the reconstructor object and the loader for a detector
c757bafd 1451
b8cd5251 1452 if (fReconstructor[iDet]) return fReconstructor[iDet];
1453
1454 // load the reconstructor object
1455 TPluginManager* pluginManager = gROOT->GetPluginManager();
1456 TString detName = fgkDetectorName[iDet];
1457 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1458 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1459
1460 if (detName == "HLT") {
1461 if (!gROOT->GetClass("AliLevel3")) {
4aa41877 1462 gSystem->Load("libAliHLTSrc.so");
1463 gSystem->Load("libAliHLTMisc.so");
1464 gSystem->Load("libAliHLTHough.so");
1465 gSystem->Load("libAliHLTComp.so");
b8cd5251 1466 }
1467 }
1468
1469 AliReconstructor* reconstructor = NULL;
1470 // first check if a plugin is defined for the reconstructor
1471 TPluginHandler* pluginHandler =
1472 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1473 // if not, add a plugin for it
1474 if (!pluginHandler) {
b8cd5251 1475 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1476 TString libs = gSystem->GetLibraries();
1477 if (libs.Contains("lib" + detName + "base.so") ||
1478 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1479 pluginManager->AddHandler("AliReconstructor", detName,
1480 recName, detName + "rec", recName + "()");
1481 } else {
1482 pluginManager->AddHandler("AliReconstructor", detName,
1483 recName, detName, recName + "()");
c757bafd 1484 }
b8cd5251 1485 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1486 }
1487 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1488 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1489 }
b8cd5251 1490 if (reconstructor) {
1491 TObject* obj = fOptions.FindObject(detName.Data());
1492 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1493 reconstructor->Init(fRunLoader);
b8cd5251 1494 fReconstructor[iDet] = reconstructor;
1495 }
1496
f08fc9f5 1497 // get or create the loader
1498 if (detName != "HLT") {
1499 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1500 if (!fLoader[iDet]) {
1501 AliConfig::Instance()
1502 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1503 detName, detName);
1504 // first check if a plugin is defined for the loader
bb0901a4 1505 pluginHandler =
f08fc9f5 1506 pluginManager->FindHandler("AliLoader", detName);
1507 // if not, add a plugin for it
1508 if (!pluginHandler) {
1509 TString loaderName = "Ali" + detName + "Loader";
1510 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1511 pluginManager->AddHandler("AliLoader", detName,
1512 loaderName, detName + "base",
1513 loaderName + "(const char*, TFolder*)");
1514 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1515 }
1516 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1517 fLoader[iDet] =
1518 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1519 fRunLoader->GetEventFolder());
1520 }
1521 if (!fLoader[iDet]) { // use default loader
1522 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1523 }
1524 if (!fLoader[iDet]) {
1525 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1526 if (fStopOnError) return NULL;
f08fc9f5 1527 } else {
1528 fRunLoader->AddLoader(fLoader[iDet]);
1529 fRunLoader->CdGAFile();
1530 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1531 fRunLoader->Write(0, TObject::kOverwrite);
1532 }
1533 }
1534 }
1535
b8cd5251 1536 return reconstructor;
c757bafd 1537}
1538
2257f27e 1539//_____________________________________________________________________________
1540Bool_t AliReconstruction::CreateVertexer()
1541{
1542// create the vertexer
1543
b8cd5251 1544 fVertexer = NULL;
1545 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1546 if (itsReconstructor) {
b8cd5251 1547 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1548 }
b8cd5251 1549 if (!fVertexer) {
815c2b38 1550 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1551 if (fStopOnError) return kFALSE;
1552 }
1553
1554 return kTRUE;
1555}
1556
24f7a148 1557//_____________________________________________________________________________
b8cd5251 1558Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1559{
f08fc9f5 1560// create the trackers
24f7a148 1561
b8cd5251 1562 TString detStr = detectors;
1563 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1564 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1565 AliReconstructor* reconstructor = GetReconstructor(iDet);
1566 if (!reconstructor) continue;
1567 TString detName = fgkDetectorName[iDet];
1f46a9ae 1568 if (detName == "HLT") {
1569 fRunHLTTracking = kTRUE;
1570 continue;
1571 }
f08fc9f5 1572
1573 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1574 if (!fTracker[iDet] && (iDet < 7)) {
1575 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1576 if (fStopOnError) return kFALSE;
1577 }
1578 }
1579
24f7a148 1580 return kTRUE;
1581}
1582
e583c30d 1583//_____________________________________________________________________________
b26c3770 1584void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1585{
1586// delete trackers and the run loader and close and delete the file
1587
b8cd5251 1588 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1589 delete fReconstructor[iDet];
1590 fReconstructor[iDet] = NULL;
1591 fLoader[iDet] = NULL;
1592 delete fTracker[iDet];
1593 fTracker[iDet] = NULL;
1594 }
1595 delete fVertexer;
1596 fVertexer = NULL;
9178838a 1597 delete fDiamondProfile;
1598 fDiamondProfile = NULL;
e583c30d 1599
1600 delete fRunLoader;
1601 fRunLoader = NULL;
b649205a 1602 delete fRawReader;
1603 fRawReader = NULL;
e583c30d 1604
1605 if (file) {
1606 file->Close();
1607 delete file;
1608 }
b26c3770 1609
1610 if (fileOld) {
1611 fileOld->Close();
1612 delete fileOld;
1613 gSystem->Unlink("AliESDs.old.root");
1614 }
e583c30d 1615}
24f7a148 1616
1617
1618//_____________________________________________________________________________
1619Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1620{
1621// read the ESD event from a file
1622
1623 if (!esd) return kFALSE;
1624 char fileName[256];
1625 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1626 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1627 if (gSystem->AccessPathName(fileName)) return kFALSE;
1628
f3a97c86 1629 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1630 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1631 TFile* file = TFile::Open(fileName);
1632 if (!file || !file->IsOpen()) {
815c2b38 1633 AliError(Form("opening %s failed", fileName));
24f7a148 1634 delete file;
1635 return kFALSE;
1636 }
1637
1638 gROOT->cd();
1639 delete esd;
1640 esd = (AliESD*) file->Get("ESD");
1641 file->Close();
1642 delete file;
1643 return kTRUE;
1644}
1645
1646//_____________________________________________________________________________
1647void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1648{
1649// write the ESD event to a file
1650
1651 if (!esd) return;
1652 char fileName[256];
1653 sprintf(fileName, "ESD_%d.%d_%s.root",
31fd97b2 1654 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
24f7a148 1655
815c2b38 1656 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1657 TFile* file = TFile::Open(fileName, "recreate");
1658 if (!file || !file->IsOpen()) {
815c2b38 1659 AliError(Form("opening %s failed", fileName));
24f7a148 1660 } else {
1661 esd->Write("ESD");
1662 file->Close();
1663 }
1664 delete file;
1665}
f3a97c86 1666
1667
1668
1669
1670//_____________________________________________________________________________
1671void AliReconstruction::CreateTag(TFile* file)
1672{
2bdb9d38 1673 /////////////
1674 //muon code//
1675 ////////////
56982dd1 1676 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1677 //Variables
56982dd1 1678 Double_t fX,fY,fZ ;
1679 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1680 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1681 Int_t fCharge;
1682 TLorentzVector fEPvector;
1683
1684 Float_t fZVertexCut = 10.0;
1685 Float_t fRhoVertexCut = 2.0;
1686
1687 Float_t fLowPtCut = 1.0;
1688 Float_t fHighPtCut = 3.0;
1689 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1690 ////////////
1691
1692 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1693
089bf903 1694 // Creates the tags for all the events in a given ESD file
4302e20f 1695 Int_t ntrack;
089bf903 1696 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1697 Int_t nPos, nNeg, nNeutr;
1698 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1699 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1700 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1701 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1702 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1703 Int_t fVertexflag;
1387d165 1704 Int_t iRunNumber = 0;
547d75a6 1705 TString fVertexName("default");
4302e20f 1706
f3a97c86 1707 AliRunTag *tag = new AliRunTag();
f3a97c86 1708 AliEventTag *evTag = new AliEventTag();
1709 TTree ttag("T","A Tree with event tags");
2bdb9d38 1710 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1711 btag->SetCompressionLevel(9);
2bdb9d38 1712
f3a97c86 1713 AliInfo(Form("Creating the tags......."));
1714
1715 if (!file || !file->IsOpen()) {
1716 AliError(Form("opening failed"));
1717 delete file;
1718 return ;
2bdb9d38 1719 }
d71aa170 1720 Int_t lastEvent = 0;
f3a97c86 1721 TTree *t = (TTree*) file->Get("esdTree");
1722 TBranch * b = t->GetBranch("ESD");
1723 AliESD *esd = 0;
1724 b->SetAddress(&esd);
1387d165 1725
d71aa170 1726 b->GetEntry(fFirstEvent);
1387d165 1727 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1728
bb0901a4 1729 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1730 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1731 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 1732 ntrack = 0;
1733 nPos = 0;
1734 nNeg = 0;
1735 nNeutr =0;
1736 nK0s = 0;
1737 nNeutrons = 0;
1738 nPi0s = 0;
1739 nGamas = 0;
1740 nProtons = 0;
1741 nKaons = 0;
1742 nPions = 0;
1743 nMuons = 0;
1744 nElectrons = 0;
1745 nCh1GeV = 0;
1746 nCh3GeV = 0;
1747 nCh10GeV = 0;
1748 nMu1GeV = 0;
1749 nMu3GeV = 0;
1750 nMu10GeV = 0;
1751 nEl1GeV = 0;
1752 nEl3GeV = 0;
1753 nEl10GeV = 0;
1754 maxPt = .0;
1755 meanPt = .0;
1756 totalP = .0;
547d75a6 1757 fVertexflag = 0;
d1a50cb5 1758
2bdb9d38 1759 b->GetEntry(iEventNumber);
1387d165 1760 iRunNumber = esd->GetRunNumber();
1761 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1762 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1763 if (!vertexIn) AliError("ESD has not defined vertex.");
1764 if (vertexIn) fVertexName = vertexIn->GetName();
1765 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1766 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1767 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1768 UInt_t status = esdTrack->GetStatus();
f3a97c86 1769
2bdb9d38 1770 //select only tracks with ITS refit
1771 if ((status&AliESDtrack::kITSrefit)==0) continue;
1772 //select only tracks with TPC refit
1773 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1774
2bdb9d38 1775 //select only tracks with the "combined PID"
1776 if ((status&AliESDtrack::kESDpid)==0) continue;
1777 Double_t p[3];
1778 esdTrack->GetPxPyPz(p);
1779 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1780 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1781 totalP += momentum;
1782 meanPt += fPt;
1783 if(fPt > maxPt) maxPt = fPt;
4302e20f 1784
2bdb9d38 1785 if(esdTrack->GetSign() > 0) {
1786 nPos++;
56982dd1 1787 if(fPt > fLowPtCut) nCh1GeV++;
1788 if(fPt > fHighPtCut) nCh3GeV++;
1789 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1790 }
1791 if(esdTrack->GetSign() < 0) {
1792 nNeg++;
56982dd1 1793 if(fPt > fLowPtCut) nCh1GeV++;
1794 if(fPt > fHighPtCut) nCh3GeV++;
1795 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1796 }
1797 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1798
2bdb9d38 1799 //PID
1800 Double_t prob[5];
1801 esdTrack->GetESDpid(prob);
4302e20f 1802
2bdb9d38 1803 Double_t rcc = 0.0;
1804 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1805 if(rcc == 0.0) continue;
1806 //Bayes' formula
1807 Double_t w[5];
1808 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1809
2bdb9d38 1810 //protons
1811 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1812 //kaons
1813 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1814 //pions
1815 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1816 //electrons
1817 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1818 nElectrons++;
56982dd1 1819 if(fPt > fLowPtCut) nEl1GeV++;
1820 if(fPt > fHighPtCut) nEl3GeV++;
1821 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1822 }
1823 ntrack++;
1824 }//track loop
1825
1826 /////////////
1827 //muon code//
1828 ////////////
1829 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1830 // loop over all reconstructed tracks (also first track of combination)
1831 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1832 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1833 if (muonTrack == 0x0) continue;
4302e20f 1834
2bdb9d38 1835 // Coordinates at vertex
56982dd1 1836 fZ = muonTrack->GetZ();
1837 fY = muonTrack->GetBendingCoor();
1838 fX = muonTrack->GetNonBendingCoor();
4302e20f 1839
56982dd1 1840 fThetaX = muonTrack->GetThetaX();
1841 fThetaY = muonTrack->GetThetaY();
4302e20f 1842
56982dd1 1843 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1844 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1845 fPxRec = fPzRec * TMath::Tan(fThetaX);
1846 fPyRec = fPzRec * TMath::Tan(fThetaY);
1847 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1848
2bdb9d38 1849 //ChiSquare of the track if needed
56982dd1 1850 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1851 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1852 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1853
2bdb9d38 1854 // total number of muons inside a vertex cut
56982dd1 1855 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1856 nMuons++;
56982dd1 1857 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1858 nMu1GeV++;
56982dd1 1859 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1860 nMu3GeV++;
56982dd1 1861 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1862 nMu10GeV++;
1863 }
1864 }
1865 }
1866 }
1867 }//muon track loop
1868
1869 // Fill the event tags
1870 if(ntrack != 0)
1871 meanPt = meanPt/ntrack;
1872
1873 evTag->SetEventId(iEventNumber+1);
547d75a6 1874 if (vertexIn) {
1875 evTag->SetVertexX(vertexIn->GetXv());
1876 evTag->SetVertexY(vertexIn->GetYv());
1877 evTag->SetVertexZ(vertexIn->GetZv());
1878 evTag->SetVertexZError(vertexIn->GetZRes());
1879 }
d1a50cb5 1880 evTag->SetVertexFlag(fVertexflag);
1881
2bdb9d38 1882 evTag->SetT0VertexZ(esd->GetT0zVertex());
1883
8bd8ac26 1884 evTag->SetTriggerMask(esd->GetTriggerMask());
1885 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1886
32a5cab4 1887 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1888 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1889 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1890 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1891 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1892 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1893
1894
1895 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1896 evTag->SetNumOfPosTracks(nPos);
1897 evTag->SetNumOfNegTracks(nNeg);
1898 evTag->SetNumOfNeutrTracks(nNeutr);
1899
1900 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1901 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1902 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1903 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1904
1905 evTag->SetNumOfProtons(nProtons);
1906 evTag->SetNumOfKaons(nKaons);
1907 evTag->SetNumOfPions(nPions);
1908 evTag->SetNumOfMuons(nMuons);
1909 evTag->SetNumOfElectrons(nElectrons);
1910 evTag->SetNumOfPhotons(nGamas);
1911 evTag->SetNumOfPi0s(nPi0s);
1912 evTag->SetNumOfNeutrons(nNeutrons);
1913 evTag->SetNumOfKaon0s(nK0s);
1914
1915 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1916 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1917 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1918 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1919 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1920 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1921 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1922 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1923 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1924
85c60a8e 1925 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1926 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1927
1928 evTag->SetTotalMomentum(totalP);
1929 evTag->SetMeanPt(meanPt);
1930 evTag->SetMaxPt(maxPt);
1931
1387d165 1932 tag->SetRunId(iInitRunNumber);
2bdb9d38 1933 tag->AddEventTag(*evTag);
1934 }
bb0901a4 1935 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 1936 else lastEvent = fLastEvent;
f3a97c86 1937
1938 ttag.Fill();
1939 tag->Clear();
1940
1941 char fileName[256];
1942 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 1943 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 1944 AliInfo(Form("writing tags to file %s", fileName));
1945 AliDebug(1, Form("writing tags to file %s", fileName));
1946
1947 TFile* ftag = TFile::Open(fileName, "recreate");
1948 ftag->cd();
1949 ttag.Write();
1950 ftag->Close();
1951 file->cd();
1952 delete tag;
f3a97c86 1953 delete evTag;
1954}
1955
a7807689 1956//_____________________________________________________________________________
1957void AliReconstruction::CreateAOD(TFile* esdFile)
1958{
1959 // do nothing for now
1960
1961 return;
1962}
1963
1964
98937d93 1965void AliReconstruction::WriteAlignmentData(AliESD* esd)
1966{
1967 // Write space-points which are then used in the alignment procedures
1968 // For the moment only ITS, TRD and TPC
1969
1970 // Load TOF clusters
d528ee75 1971 if (fTracker[3]){
1972 fLoader[3]->LoadRecPoints("read");
1973 TTree* tree = fLoader[3]->TreeR();
1974 if (!tree) {
1975 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1976 return;
1977 }
1978 fTracker[3]->LoadClusters(tree);
98937d93 1979 }
98937d93 1980 Int_t ntracks = esd->GetNumberOfTracks();
1981 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1982 {
1983 AliESDtrack *track = esd->GetTrack(itrack);
1984 Int_t nsp = 0;
ef7253ac 1985 Int_t idx[200];
98937d93 1986 for (Int_t iDet = 3; iDet >= 0; iDet--)
1987 nsp += track->GetNcls(iDet);
1988 if (nsp) {
1989 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1990 track->SetTrackPointArray(sp);
1991 Int_t isptrack = 0;
1992 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1993 AliTracker *tracker = fTracker[iDet];
1994 if (!tracker) continue;
1995 Int_t nspdet = track->GetNcls(iDet);
98937d93 1996 if (nspdet <= 0) continue;
1997 track->GetClusters(iDet,idx);
1998 AliTrackPoint p;
1999 Int_t isp = 0;
2000 Int_t isp2 = 0;
2001 while (isp < nspdet) {
2002 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 2003 const Int_t kNTPCmax = 159;
2004 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 2005 if (!isvalid) continue;
2006 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2007 }
98937d93 2008 }
2009 }
2010 }
d528ee75 2011 if (fTracker[3]){
2012 fTracker[3]->UnloadClusters();
2013 fLoader[3]->UnloadRecPoints();
2014 }
98937d93 2015}