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