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