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