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