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