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