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