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