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