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