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