Corrected position of the compensation dipole.
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
CommitLineData
596a855f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// class for running the reconstruction //
21// //
22// Clusters and tracks are created for all detectors and all events by //
23// typing: //
24// //
25// AliReconstruction rec; //
26// rec.Run(); //
27// //
28// The Run method returns kTRUE in case of successful execution. //
c71de921 29// //
30// If the input to the reconstruction are not simulated digits but raw data, //
31// this can be specified by an argument of the Run method or by the method //
32// //
33// rec.SetInput("..."); //
34// //
35// The input formats and the corresponding argument are: //
36// - DDL raw data files: directory name, ends with "/" //
37// - raw data root file: root file name, extension ".root" //
38// - raw data DATE file: DATE file name, any other non-empty string //
39// - MC root files : empty string, default //
40// //
b26c3770 41// By default all events are reconstructed. The reconstruction can be //
42// limited to a range of events by giving the index of the first and the //
43// last event as an argument to the Run method or by calling //
44// //
45// rec.SetEventRange(..., ...); //
46// //
47// The index -1 (default) can be used for the last event to indicate no //
48// upper limit of the event range. //
49// //
596a855f 50// The name of the galice file can be changed from the default //
e583c30d 51// "galice.root" by passing it as argument to the AliReconstruction //
52// constructor or by //
596a855f 53// //
54// rec.SetGAliceFile("..."); //
55// //
59697224 56// The local reconstruction can be switched on or off for individual //
57// detectors by //
596a855f 58// //
59697224 59// rec.SetRunLocalReconstruction("..."); //
596a855f 60// //
61// The argument is a (case sensitive) string with the names of the //
62// detectors separated by a space. The special string "ALL" selects all //
63// available detectors. This is the default. //
64// //
c71de921 65// The reconstruction of the primary vertex position can be switched off by //
66// //
67// rec.SetRunVertexFinder(kFALSE); //
68// //
b8cd5251 69// The tracking and the creation of ESD tracks can be switched on for //
70// selected detectors by //
596a855f 71// //
b8cd5251 72// rec.SetRunTracking("..."); //
596a855f 73// //
c84a5e9e 74// Uniform/nonuniform field tracking switches (default: uniform field) //
75// //
1d99986f 76// rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
c84a5e9e 77// //
596a855f 78// The filling of additional ESD information can be steered by //
79// //
80// rec.SetFillESD("..."); //
81// //
b8cd5251 82// Again, for both methods the string specifies the list of detectors. //
83// The default is "ALL". //
84// //
85// The call of the shortcut method //
86// //
87// rec.SetRunReconstruction("..."); //
88// //
89// is equivalent to calling SetRunLocalReconstruction, SetRunTracking and //
90// SetFillESD with the same detector selecting string as argument. //
596a855f 91// //
c71de921 92// The reconstruction requires digits or raw data as input. For the creation //
93// of digits and raw data have a look at the class AliSimulation. //
596a855f 94// //
24f7a148 95// For debug purposes the method SetCheckPointLevel can be used. If the //
96// argument is greater than 0, files with ESD events will be written after //
97// selected steps of the reconstruction for each event: //
98// level 1: after tracking and after filling of ESD (final) //
99// level 2: in addition after each tracking step //
100// level 3: in addition after the filling of ESD for each detector //
101// If a final check point file exists for an event, this event will be //
102// skipped in the reconstruction. The tracking and the filling of ESD for //
103// a detector will be skipped as well, if the corresponding check point //
104// file exists. The ESD event will then be loaded from the file instead. //
105// //
596a855f 106///////////////////////////////////////////////////////////////////////////////
107
024a7e64 108#include <TArrayF.h>
109#include <TFile.h>
110#include <TSystem.h>
111#include <TROOT.h>
112#include <TPluginManager.h>
fd46e2d2 113#include <TStopwatch.h>
3103d196 114#include <TGeoManager.h>
2bdb9d38 115#include <TLorentzVector.h>
596a855f 116
117#include "AliReconstruction.h"
b8cd5251 118#include "AliReconstructor.h"
815c2b38 119#include "AliLog.h"
596a855f 120#include "AliRunLoader.h"
121#include "AliRun.h"
b649205a 122#include "AliRawReaderFile.h"
123#include "AliRawReaderDate.h"
124#include "AliRawReaderRoot.h"
001397cd 125#include "AliRawEventHeaderBase.h"
596a855f 126#include "AliESD.h"
1d99986f 127#include "AliESDfriend.h"
2257f27e 128#include "AliESDVertex.h"
32e449be 129#include "AliMultiplicity.h"
c84a5e9e 130#include "AliTracker.h"
2257f27e 131#include "AliVertexer.h"
c5e3e5d1 132#include "AliVertexerTracks.h"
596a855f 133#include "AliHeader.h"
134#include "AliGenEventHeader.h"
b26c3770 135#include "AliPID.h"
596a855f 136#include "AliESDpid.h"
ff8bb5ae 137#include "AliESDtrack.h"
f3a97c86 138
139#include "AliRunTag.h"
f3a97c86 140#include "AliDetectorTag.h"
141#include "AliEventTag.h"
142
98937d93 143#include "AliTrackPointArray.h"
b0314964 144#include "AliCDBManager.h"
6bae477a 145#include "AliCDBEntry.h"
146#include "AliAlignObj.h"
f3a97c86 147
b647652d 148#include "AliCentralTrigger.h"
149#include "AliCTPRawStream.h"
150
596a855f 151ClassImp(AliReconstruction)
152
153
154//_____________________________________________________________________________
482070f2 155const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
c757bafd 156
157//_____________________________________________________________________________
024cf675 158AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 159 const char* name, const char* title) :
160 TNamed(name, title),
161
c84a5e9e 162 fUniformField(kTRUE),
2257f27e 163 fRunVertexFinder(kTRUE),
1f46a9ae 164 fRunHLTTracking(kFALSE),
1d99986f 165 fStopOnError(kFALSE),
166 fWriteAlignmentData(kFALSE),
167 fWriteESDfriend(kFALSE),
b647652d 168 fFillTriggerESD(kTRUE),
1d99986f 169
170 fRunLocalReconstruction("ALL"),
b8cd5251 171 fRunTracking("ALL"),
e583c30d 172 fFillESD("ALL"),
173 fGAliceFileName(gAliceFilename),
b649205a 174 fInput(""),
35042093 175 fEquipIdMap(""),
b26c3770 176 fFirstEvent(0),
177 fLastEvent(-1),
24f7a148 178 fCheckPointLevel(0),
b8cd5251 179 fOptions(),
6bae477a 180 fLoadAlignFromCDB(kTRUE),
181 fLoadAlignData("ALL"),
e583c30d 182
183 fRunLoader(NULL),
b649205a 184 fRawReader(NULL),
b8cd5251 185
98937d93 186 fVertexer(NULL),
9178838a 187 fDiamondProfile(NULL),
98937d93 188
6bae477a 189 fAlignObjArray(NULL),
ec92bee0 190 fCDBUri(cdbUri),
191 fSpecCDBUri()
596a855f 192{
193// create reconstruction object with default parameters
b8cd5251 194
195 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
196 fReconstructor[iDet] = NULL;
197 fLoader[iDet] = NULL;
198 fTracker[iDet] = NULL;
199 }
e47c4c2e 200 AliPID pid;
596a855f 201}
202
203//_____________________________________________________________________________
204AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 205 TNamed(rec),
206
c84a5e9e 207 fUniformField(rec.fUniformField),
2257f27e 208 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 209 fRunHLTTracking(rec.fRunHLTTracking),
1d99986f 210 fStopOnError(rec.fStopOnError),
211 fWriteAlignmentData(rec.fWriteAlignmentData),
212 fWriteESDfriend(rec.fWriteESDfriend),
b647652d 213 fFillTriggerESD(rec.fFillTriggerESD),
1d99986f 214
215 fRunLocalReconstruction(rec.fRunLocalReconstruction),
e583c30d 216 fRunTracking(rec.fRunTracking),
217 fFillESD(rec.fFillESD),
218 fGAliceFileName(rec.fGAliceFileName),
b649205a 219 fInput(rec.fInput),
35042093 220 fEquipIdMap(rec.fEquipIdMap),
b26c3770 221 fFirstEvent(rec.fFirstEvent),
222 fLastEvent(rec.fLastEvent),
24f7a148 223 fCheckPointLevel(0),
b8cd5251 224 fOptions(),
6bae477a 225 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
226 fLoadAlignData(rec.fLoadAlignData),
e583c30d 227
228 fRunLoader(NULL),
b649205a 229 fRawReader(NULL),
b8cd5251 230
98937d93 231 fVertexer(NULL),
9178838a 232 fDiamondProfile(NULL),
98937d93 233
6bae477a 234 fAlignObjArray(rec.fAlignObjArray),
ec92bee0 235 fCDBUri(rec.fCDBUri),
236 fSpecCDBUri()
596a855f 237{
238// copy constructor
239
ec92bee0 240 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
efd2085e 241 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
242 }
b8cd5251 243 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
244 fReconstructor[iDet] = NULL;
245 fLoader[iDet] = NULL;
246 fTracker[iDet] = NULL;
247 }
ec92bee0 248 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
249 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
250 }
596a855f 251}
252
253//_____________________________________________________________________________
254AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
255{
256// assignment operator
257
258 this->~AliReconstruction();
259 new(this) AliReconstruction(rec);
260 return *this;
261}
262
263//_____________________________________________________________________________
264AliReconstruction::~AliReconstruction()
265{
266// clean up
267
e583c30d 268 CleanUp();
efd2085e 269 fOptions.Delete();
ec92bee0 270 fSpecCDBUri.Delete();
596a855f 271}
272
024cf675 273//_____________________________________________________________________________
274void AliReconstruction::InitCDBStorage()
275{
276// activate a default CDB storage
277// First check if we have any CDB storage set, because it is used
278// to retrieve the calibration and alignment constants
279
280 AliCDBManager* man = AliCDBManager::Instance();
ec92bee0 281 if (man->IsDefaultStorageSet())
024cf675 282 {
283 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 284 AliWarning("Default CDB storage has been already set !");
285 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
024cf675 286 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 287 fCDBUri = "";
288 }
289 else {
b8ec52f6 290 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
291 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
292 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 293 man->SetDefaultStorage(fCDBUri);
294 }
295
296 // Now activate the detector specific CDB storage locations
c3a7b59a 297 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
298 TObject* obj = fSpecCDBUri[i];
299 if (!obj) continue;
b8ec52f6 300 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
301 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
302 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
c3a7b59a 303 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
ec92bee0 304 }
727d0766 305 man->Print();
024cf675 306}
307
308//_____________________________________________________________________________
309void AliReconstruction::SetDefaultStorage(const char* uri) {
ec92bee0 310// Store the desired default CDB storage location
311// Activate it later within the Run() method
024cf675 312
ec92bee0 313 fCDBUri = uri;
024cf675 314
315}
316
317//_____________________________________________________________________________
c3a7b59a 318void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
ec92bee0 319// Store a detector-specific CDB storage location
320// Activate it later within the Run() method
024cf675 321
c3a7b59a 322 AliCDBPath aPath(calibType);
323 if(!aPath.IsValid()){
324 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
325 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
326 if(!strcmp(calibType, fgkDetectorName[iDet])) {
327 aPath.SetPath(Form("%s/*", calibType));
328 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
329 break;
330 }
331 }
332 if(!aPath.IsValid()){
333 AliError(Form("Not a valid path or detector: %s", calibType));
334 return;
335 }
336 }
337
338 // check that calibType refers to a "valid" detector name
339 Bool_t isDetector = kFALSE;
340 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
341 TString detName = fgkDetectorName[iDet];
342 if(aPath.GetLevel0() == detName) {
343 isDetector = kTRUE;
344 break;
345 }
346 }
347
348 if(!isDetector) {
349 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
350 return;
351 }
352
353 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
ec92bee0 354 if (obj) fSpecCDBUri.Remove(obj);
c3a7b59a 355 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
024cf675 356
357}
358
6bae477a 359//_____________________________________________________________________________
360Bool_t AliReconstruction::SetRunNumber()
361{
362 // The method is called in Run() in order
363 // to set a correct run number.
364 // In case of raw data reconstruction the
365 // run number is taken from the raw data header
366
367 if(AliCDBManager::Instance()->GetRun() < 0) {
368 if (!fRunLoader) {
369 AliError("No run loader is found !");
370 return kFALSE;
371 }
372 // read run number from gAlice
ec92bee0 373 if(fRunLoader->GetAliRun())
374 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
375 else {
376 if(fRawReader) {
377 if(fRawReader->NextEvent()) {
378 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
379 fRawReader->RewindEvents();
380 }
381 else {
382 AliError("No raw-data events found !");
383 return kFALSE;
384 }
385 }
386 else {
387 AliError("Neither gAlice nor RawReader objects are found !");
388 return kFALSE;
389 }
390 }
391 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
6bae477a 392 }
393 return kTRUE;
394}
395
396//_____________________________________________________________________________
397Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
398{
399 // Read collection of alignment objects (AliAlignObj derived) saved
400 // in the TClonesArray ClArrayName and apply them to the geometry
401 // manager singleton.
402 //
403 alObjArray->Sort();
404 Int_t nvols = alObjArray->GetEntriesFast();
405
dc0984f8 406 Bool_t flag = kTRUE;
407
6bae477a 408 for(Int_t j=0; j<nvols; j++)
409 {
410 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
dc0984f8 411 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
6bae477a 412 }
413
414 if (AliDebugLevelClass() >= 1) {
e0625db4 415 gGeoManager->GetTopNode()->CheckOverlaps(20);
6bae477a 416 TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
417 if(ovexlist->GetEntriesFast()){
418 AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
419 }
420 }
421
dc0984f8 422 return flag;
6bae477a 423
424}
425
426//_____________________________________________________________________________
427Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
428{
429 // Fills array of single detector's alignable objects from CDB
430
431 AliDebug(2, Form("Loading alignment data for detector: %s",detName));
432
433 AliCDBEntry *entry;
434
435 AliCDBPath path(detName,"Align","Data");
436
437 entry=AliCDBManager::Instance()->Get(path.GetPath());
438 if(!entry){
439 AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
440 return kFALSE;
441 }
442 entry->SetOwner(1);
443 TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
444 alignArray->SetOwner(0);
445 AliDebug(2,Form("Found %d alignment objects for %s",
446 alignArray->GetEntries(),detName));
447
448 AliAlignObj *alignObj=0;
449 TIter iter(alignArray);
450
451 // loop over align objects in detector
452 while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
453 fAlignObjArray->Add(alignObj);
454 }
455 // delete entry --- Don't delete, it is cached!
456
457 AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
458 return kTRUE;
459
460}
461
462//_____________________________________________________________________________
463Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
464{
465 // Read the alignment objects from CDB.
466 // Each detector is supposed to have the
467 // alignment objects in DET/Align/Data CDB path.
468 // All the detector objects are then collected,
469 // sorted by geometry level (starting from ALIC) and
470 // then applied to the TGeo geometry.
471 // Finally an overlaps check is performed.
472
473 // Load alignment data from CDB and fill fAlignObjArray
474 if(fLoadAlignFromCDB){
475 if(!fAlignObjArray) fAlignObjArray = new TObjArray();
476
477 //fAlignObjArray->RemoveAll();
478 fAlignObjArray->Clear();
479 fAlignObjArray->SetOwner(0);
480
481 TString detStr = detectors;
482 TString dataNotLoaded="";
483 TString dataLoaded="";
484
485 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
486 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
487 if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
488 dataNotLoaded += fgkDetectorName[iDet];
489 dataNotLoaded += " ";
490 } else {
491 dataLoaded += fgkDetectorName[iDet];
492 dataLoaded += " ";
493 }
494 } // end loop over detectors
495
496 if ((detStr.CompareTo("ALL") == 0)) detStr = "";
497 dataNotLoaded += detStr;
7250c343 498 if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
6bae477a 499 dataLoaded.Data()));
7250c343 500 if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
6bae477a 501 dataNotLoaded.Data()));
502 } // fLoadAlignFromCDB flag
503
504 // Check if the array with alignment objects was
505 // provided by the user. If yes, apply the objects
506 // to the present TGeo geometry
507 if (fAlignObjArray) {
508 if (gGeoManager && gGeoManager->IsClosed()) {
509 if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
dc0984f8 510 AliError("The misalignment of one or more volumes failed!"
511 "Compare the list of simulated detectors and the list of detector alignment data!");
6bae477a 512 return kFALSE;
513 }
514 }
515 else {
516 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
517 return kFALSE;
518 }
519 }
520
521 return kTRUE;
522}
596a855f 523
524//_____________________________________________________________________________
525void AliReconstruction::SetGAliceFile(const char* fileName)
526{
527// set the name of the galice file
528
529 fGAliceFileName = fileName;
530}
531
efd2085e 532//_____________________________________________________________________________
533void AliReconstruction::SetOption(const char* detector, const char* option)
534{
535// set options for the reconstruction of a detector
536
537 TObject* obj = fOptions.FindObject(detector);
538 if (obj) fOptions.Remove(obj);
539 fOptions.Add(new TNamed(detector, option));
540}
541
596a855f 542
543//_____________________________________________________________________________
b26c3770 544Bool_t AliReconstruction::Run(const char* input,
545 Int_t firstEvent, Int_t lastEvent)
596a855f 546{
547// run the reconstruction
548
b649205a 549 // set the input
550 if (!input) input = fInput.Data();
551 TString fileName(input);
552 if (fileName.EndsWith("/")) {
553 fRawReader = new AliRawReaderFile(fileName);
554 } else if (fileName.EndsWith(".root")) {
555 fRawReader = new AliRawReaderRoot(fileName);
556 } else if (!fileName.IsNull()) {
557 fRawReader = new AliRawReaderDate(fileName);
558 fRawReader->SelectEvents(7);
559 }
35042093 560 if (!fEquipIdMap.IsNull() && fRawReader)
561 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
562
b649205a 563
f08fc9f5 564 // get the run loader
565 if (!InitRunLoader()) return kFALSE;
596a855f 566
ec92bee0 567 // Initialize the CDB storage
568 InitCDBStorage();
569
6bae477a 570 // Set run number in CDBManager (if it is not already set by the user)
571 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
572
573 // Import ideal TGeo geometry and apply misalignment
574 if (!gGeoManager) {
575 TString geom(gSystem->DirName(fGAliceFileName));
576 geom += "/geometry.root";
577 TGeoManager::Import(geom.Data());
578 if (!gGeoManager) if (fStopOnError) return kFALSE;
579 }
580 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
581
596a855f 582 // local reconstruction
59697224 583 if (!fRunLocalReconstruction.IsNull()) {
584 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 585 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 586 }
587 }
b26c3770 588// if (!fRunVertexFinder && fRunTracking.IsNull() &&
589// fFillESD.IsNull()) return kTRUE;
2257f27e 590
591 // get vertexer
592 if (fRunVertexFinder && !CreateVertexer()) {
593 if (fStopOnError) {
594 CleanUp();
595 return kFALSE;
596 }
597 }
596a855f 598
f08fc9f5 599 // get trackers
b8cd5251 600 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 601 if (fStopOnError) {
602 CleanUp();
603 return kFALSE;
604 }
596a855f 605 }
24f7a148 606
9db3a215 607
608 TStopwatch stopwatch;
609 stopwatch.Start();
610
b26c3770 611 // get the possibly already existing ESD file and tree
1f46a9ae 612 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
b26c3770 613 TFile* fileOld = NULL;
1f46a9ae 614 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 615 if (!gSystem->AccessPathName("AliESDs.root")){
616 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
617 fileOld = TFile::Open("AliESDs.old.root");
618 if (fileOld && fileOld->IsOpen()) {
619 treeOld = (TTree*) fileOld->Get("esdTree");
620 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
1f46a9ae 621 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
622 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
b26c3770 623 }
624 }
625
36711aa4 626 // create the ESD output file and tree
596a855f 627 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
628 if (!file->IsOpen()) {
815c2b38 629 AliError("opening AliESDs.root failed");
b26c3770 630 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 631 }
36711aa4 632 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
633 tree->Branch("ESD", "AliESD", &esd);
1f46a9ae 634 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
635 hlttree->Branch("ESD", "AliESD", &hltesd);
636 delete esd; delete hltesd;
637 esd = NULL; hltesd = NULL;
1d99986f 638
500d54ab 639 // create the branch with ESD additions
640 AliESDfriend *esdf=0;
1d99986f 641 if (fWriteESDfriend) {
500d54ab 642 TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
643 br->SetFile("AliESDfriends.root");
1d99986f 644 }
645
c5e3e5d1 646 AliVertexerTracks tVertexer;
9178838a 647 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
c5e3e5d1 648
596a855f 649 // loop over events
b649205a 650 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 651
596a855f 652 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 653 if (fRawReader) fRawReader->NextEvent();
654 if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
655 // copy old ESD to the new one
656 if (treeOld) {
657 treeOld->SetBranchAddress("ESD", &esd);
658 treeOld->GetEntry(iEvent);
659 }
660 tree->Fill();
1f46a9ae 661 if (hlttreeOld) {
662 hlttreeOld->SetBranchAddress("ESD", &hltesd);
663 hlttreeOld->GetEntry(iEvent);
664 }
665 hlttree->Fill();
b26c3770 666 continue;
667 }
668
815c2b38 669 AliInfo(Form("processing event %d", iEvent));
596a855f 670 fRunLoader->GetEvent(iEvent);
24f7a148 671
672 char fileName[256];
673 sprintf(fileName, "ESD_%d.%d_final.root",
f08fc9f5 674 fRunLoader->GetHeader()->GetRun(),
675 fRunLoader->GetHeader()->GetEventNrInRun());
24f7a148 676 if (!gSystem->AccessPathName(fileName)) continue;
677
b26c3770 678 // local reconstruction
679 if (!fRunLocalReconstruction.IsNull()) {
680 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
681 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
682 }
683 }
684
1f46a9ae 685 esd = new AliESD; hltesd = new AliESD;
f08fc9f5 686 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 687 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
f08fc9f5 688 esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
1f46a9ae 689 hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
d6ee376f 690
691 // Set magnetic field from the tracker
692 esd->SetMagneticField(AliTracker::GetBz());
693 hltesd->SetMagneticField(AliTracker::GetBz());
596a855f 694
2257f27e 695 // vertex finder
696 if (fRunVertexFinder) {
697 if (!ReadESD(esd, "vertex")) {
698 if (!RunVertexFinder(esd)) {
b26c3770 699 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 700 }
701 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
702 }
703 }
704
1f46a9ae 705 // HLT tracking
706 if (!fRunTracking.IsNull()) {
707 if (fRunHLTTracking) {
708 hltesd->SetVertex(esd->GetVertex());
709 if (!RunHLTTracking(hltesd)) {
710 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
711 }
712 }
713 }
714
596a855f 715 // barrel tracking
b8cd5251 716 if (!fRunTracking.IsNull()) {
24f7a148 717 if (!ReadESD(esd, "tracking")) {
718 if (!RunTracking(esd)) {
b26c3770 719 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 720 }
721 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
596a855f 722 }
723 }
724
725 // fill ESD
726 if (!fFillESD.IsNull()) {
727 if (!FillESD(esd, fFillESD)) {
b26c3770 728 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 729 }
730 }
001397cd 731 // fill Event header information from the RawEventHeader
732 if (fRawReader){FillRawEventHeaderESD(esd);}
596a855f 733
734 // combined PID
735 AliESDpid::MakePID(esd);
24f7a148 736 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 737
b647652d 738 if (fFillTriggerESD) {
739 if (!ReadESD(esd, "trigger")) {
740 if (!FillTriggerESD(esd)) {
741 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
742 }
743 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
744 }
745 }
746
a28195f7 747 esd->SetPrimaryVertex(tVertexer.FindPrimaryVertex(esd));
c5e3e5d1 748
596a855f 749 // write ESD
1d99986f 750 if (fWriteESDfriend) {
d75007f6 751 esdf=new AliESDfriend();
752 esd->GetESDfriend(esdf);
1d99986f 753 }
500d54ab 754 tree->Fill();
755
756 // write HLT ESD
757 hlttree->Fill();
1d99986f 758
f3a97c86 759 if (fCheckPointLevel > 0) WriteESD(esd, "final");
760
500d54ab 761 delete esd; delete esdf; delete hltesd;
762 esd = NULL; esdf=NULL; hltesd = NULL;
596a855f 763 }
764
9db3a215 765 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
766 stopwatch.RealTime(),stopwatch.CpuTime()));
767
36711aa4 768 file->cd();
500d54ab 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
823//_____________________________________________________________________________
b26c3770 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
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
957//_____________________________________________________________________________
1f46a9ae 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
1017//_____________________________________________________________________________
2257f27e 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
1306//_____________________________________________________________________________
f08fc9f5 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
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
1470//_____________________________________________________________________________
2257f27e 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
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
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 }
1651 Int_t firstEvent = 0,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
1657 b->GetEntry(0);
1658 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1659
089bf903 1660 Int_t iNumberOfEvents = b->GetEntries();
2bdb9d38 1661 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
1662 ntrack = 0;
1663 nPos = 0;
1664 nNeg = 0;
1665 nNeutr =0;
1666 nK0s = 0;
1667 nNeutrons = 0;
1668 nPi0s = 0;
1669 nGamas = 0;
1670 nProtons = 0;
1671 nKaons = 0;
1672 nPions = 0;
1673 nMuons = 0;
1674 nElectrons = 0;
1675 nCh1GeV = 0;
1676 nCh3GeV = 0;
1677 nCh10GeV = 0;
1678 nMu1GeV = 0;
1679 nMu3GeV = 0;
1680 nMu10GeV = 0;
1681 nEl1GeV = 0;
1682 nEl3GeV = 0;
1683 nEl10GeV = 0;
1684 maxPt = .0;
1685 meanPt = .0;
1686 totalP = .0;
547d75a6 1687 fVertexflag = 0;
d1a50cb5 1688
2bdb9d38 1689 b->GetEntry(iEventNumber);
1387d165 1690 iRunNumber = esd->GetRunNumber();
1691 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1692 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1693 if (!vertexIn) AliError("ESD has not defined vertex.");
1694 if (vertexIn) fVertexName = vertexIn->GetName();
1695 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1696 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1697 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1698 UInt_t status = esdTrack->GetStatus();
f3a97c86 1699
2bdb9d38 1700 //select only tracks with ITS refit
1701 if ((status&AliESDtrack::kITSrefit)==0) continue;
1702 //select only tracks with TPC refit
1703 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1704
2bdb9d38 1705 //select only tracks with the "combined PID"
1706 if ((status&AliESDtrack::kESDpid)==0) continue;
1707 Double_t p[3];
1708 esdTrack->GetPxPyPz(p);
1709 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1710 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1711 totalP += momentum;
1712 meanPt += fPt;
1713 if(fPt > maxPt) maxPt = fPt;
4302e20f 1714
2bdb9d38 1715 if(esdTrack->GetSign() > 0) {
1716 nPos++;
56982dd1 1717 if(fPt > fLowPtCut) nCh1GeV++;
1718 if(fPt > fHighPtCut) nCh3GeV++;
1719 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1720 }
1721 if(esdTrack->GetSign() < 0) {
1722 nNeg++;
56982dd1 1723 if(fPt > fLowPtCut) nCh1GeV++;
1724 if(fPt > fHighPtCut) nCh3GeV++;
1725 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1726 }
1727 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1728
2bdb9d38 1729 //PID
1730 Double_t prob[5];
1731 esdTrack->GetESDpid(prob);
4302e20f 1732
2bdb9d38 1733 Double_t rcc = 0.0;
1734 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1735 if(rcc == 0.0) continue;
1736 //Bayes' formula
1737 Double_t w[5];
1738 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1739
2bdb9d38 1740 //protons
1741 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1742 //kaons
1743 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1744 //pions
1745 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1746 //electrons
1747 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1748 nElectrons++;
56982dd1 1749 if(fPt > fLowPtCut) nEl1GeV++;
1750 if(fPt > fHighPtCut) nEl3GeV++;
1751 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1752 }
1753 ntrack++;
1754 }//track loop
1755
1756 /////////////
1757 //muon code//
1758 ////////////
1759 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1760 // loop over all reconstructed tracks (also first track of combination)
1761 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1762 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1763 if (muonTrack == 0x0) continue;
4302e20f 1764
2bdb9d38 1765 // Coordinates at vertex
56982dd1 1766 fZ = muonTrack->GetZ();
1767 fY = muonTrack->GetBendingCoor();
1768 fX = muonTrack->GetNonBendingCoor();
4302e20f 1769
56982dd1 1770 fThetaX = muonTrack->GetThetaX();
1771 fThetaY = muonTrack->GetThetaY();
4302e20f 1772
56982dd1 1773 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1774 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1775 fPxRec = fPzRec * TMath::Tan(fThetaX);
1776 fPyRec = fPzRec * TMath::Tan(fThetaY);
1777 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1778
2bdb9d38 1779 //ChiSquare of the track if needed
56982dd1 1780 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1781 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1782 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1783
2bdb9d38 1784 // total number of muons inside a vertex cut
56982dd1 1785 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1786 nMuons++;
56982dd1 1787 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1788 nMu1GeV++;
56982dd1 1789 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1790 nMu3GeV++;
56982dd1 1791 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1792 nMu10GeV++;
1793 }
1794 }
1795 }
1796 }
1797 }//muon track loop
1798
1799 // Fill the event tags
1800 if(ntrack != 0)
1801 meanPt = meanPt/ntrack;
1802
1803 evTag->SetEventId(iEventNumber+1);
547d75a6 1804 if (vertexIn) {
1805 evTag->SetVertexX(vertexIn->GetXv());
1806 evTag->SetVertexY(vertexIn->GetYv());
1807 evTag->SetVertexZ(vertexIn->GetZv());
1808 evTag->SetVertexZError(vertexIn->GetZRes());
1809 }
d1a50cb5 1810 evTag->SetVertexFlag(fVertexflag);
1811
2bdb9d38 1812 evTag->SetT0VertexZ(esd->GetT0zVertex());
1813
8bd8ac26 1814 evTag->SetTriggerMask(esd->GetTriggerMask());
1815 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1816
32a5cab4 1817 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1818 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1819 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1820 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1821 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1822 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1823
1824
1825 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1826 evTag->SetNumOfPosTracks(nPos);
1827 evTag->SetNumOfNegTracks(nNeg);
1828 evTag->SetNumOfNeutrTracks(nNeutr);
1829
1830 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1831 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1832 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1833 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1834
1835 evTag->SetNumOfProtons(nProtons);
1836 evTag->SetNumOfKaons(nKaons);
1837 evTag->SetNumOfPions(nPions);
1838 evTag->SetNumOfMuons(nMuons);
1839 evTag->SetNumOfElectrons(nElectrons);
1840 evTag->SetNumOfPhotons(nGamas);
1841 evTag->SetNumOfPi0s(nPi0s);
1842 evTag->SetNumOfNeutrons(nNeutrons);
1843 evTag->SetNumOfKaon0s(nK0s);
1844
1845 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1846 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1847 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1848 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1849 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1850 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1851 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1852 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1853 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1854
85c60a8e 1855 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1856 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1857
1858 evTag->SetTotalMomentum(totalP);
1859 evTag->SetMeanPt(meanPt);
1860 evTag->SetMaxPt(maxPt);
1861
1387d165 1862 tag->SetRunId(iInitRunNumber);
2bdb9d38 1863 tag->AddEventTag(*evTag);
1864 }
089bf903 1865 lastEvent = iNumberOfEvents;
f3a97c86 1866
1867 ttag.Fill();
1868 tag->Clear();
1869
1870 char fileName[256];
1871 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1872 tag->GetRunId(),firstEvent,lastEvent );
1873 AliInfo(Form("writing tags to file %s", fileName));
1874 AliDebug(1, Form("writing tags to file %s", fileName));
1875
1876 TFile* ftag = TFile::Open(fileName, "recreate");
1877 ftag->cd();
1878 ttag.Write();
1879 ftag->Close();
1880 file->cd();
1881 delete tag;
f3a97c86 1882 delete evTag;
1883}
1884
98937d93 1885void AliReconstruction::WriteAlignmentData(AliESD* esd)
1886{
1887 // Write space-points which are then used in the alignment procedures
1888 // For the moment only ITS, TRD and TPC
1889
1890 // Load TOF clusters
d528ee75 1891 if (fTracker[3]){
1892 fLoader[3]->LoadRecPoints("read");
1893 TTree* tree = fLoader[3]->TreeR();
1894 if (!tree) {
1895 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1896 return;
1897 }
1898 fTracker[3]->LoadClusters(tree);
98937d93 1899 }
98937d93 1900 Int_t ntracks = esd->GetNumberOfTracks();
1901 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1902 {
1903 AliESDtrack *track = esd->GetTrack(itrack);
1904 Int_t nsp = 0;
ef7253ac 1905 Int_t idx[200];
98937d93 1906 for (Int_t iDet = 3; iDet >= 0; iDet--)
1907 nsp += track->GetNcls(iDet);
1908 if (nsp) {
1909 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1910 track->SetTrackPointArray(sp);
1911 Int_t isptrack = 0;
1912 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1913 AliTracker *tracker = fTracker[iDet];
1914 if (!tracker) continue;
1915 Int_t nspdet = track->GetNcls(iDet);
98937d93 1916 if (nspdet <= 0) continue;
1917 track->GetClusters(iDet,idx);
1918 AliTrackPoint p;
1919 Int_t isp = 0;
1920 Int_t isp2 = 0;
1921 while (isp < nspdet) {
1922 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 1923 const Int_t kNTPCmax = 159;
1924 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 1925 if (!isvalid) continue;
1926 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1927 }
98937d93 1928 }
1929 }
1930 }
d528ee75 1931 if (fTracker[3]){
1932 fTracker[3]->UnloadClusters();
1933 fLoader[3]->UnloadRecPoints();
1934 }
98937d93 1935}