]> 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
29-jun-2007 NvE Unused variable dist2 removed from IceDwalk::Exec to prevent
[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 "AliGeomManager.h"
155#include "AliTrackPointArray.h"
156#include "AliCDBManager.h"
157#include "AliCDBEntry.h"
158#include "AliAlignObj.h"
159
160#include "AliCentralTrigger.h"
161#include "AliCTPRawStream.h"
162
163#include "AliAODEvent.h"
164#include "AliAODHeader.h"
165#include "AliAODTrack.h"
166#include "AliAODVertex.h"
167#include "AliAODCluster.h"
168
169
170ClassImp(AliReconstruction)
171
172
173//_____________________________________________________________________________
174const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
175
176//_____________________________________________________________________________
177AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
178 const char* name, const char* title) :
179 TNamed(name, title),
180
181 fUniformField(kTRUE),
182 fRunVertexFinder(kTRUE),
183 fRunHLTTracking(kFALSE),
184 fRunMuonTracking(kFALSE),
185 fStopOnError(kFALSE),
186 fWriteAlignmentData(kFALSE),
187 fWriteESDfriend(kFALSE),
188 fWriteAOD(kFALSE),
189 fFillTriggerESD(kTRUE),
190
191 fRunLocalReconstruction("ALL"),
192 fRunTracking("ALL"),
193 fFillESD("ALL"),
194 fGAliceFileName(gAliceFilename),
195 fInput(""),
196 fEquipIdMap(""),
197 fFirstEvent(0),
198 fLastEvent(-1),
199 fNumberOfEventsPerFile(1),
200 fCheckPointLevel(0),
201 fOptions(),
202 fLoadAlignFromCDB(kTRUE),
203 fLoadAlignData("ALL"),
204 fESDPar(""),
205
206 fRunLoader(NULL),
207 fRawReader(NULL),
208
209 fVertexer(NULL),
210 fDiamondProfile(NULL),
211
212 fAlignObjArray(NULL),
213 fCDBUri(cdbUri),
214 fSpecCDBUri()
215{
216// create reconstruction object with default parameters
217
218 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
219 fReconstructor[iDet] = NULL;
220 fLoader[iDet] = NULL;
221 fTracker[iDet] = NULL;
222 }
223 AliPID pid;
224}
225
226//_____________________________________________________________________________
227AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
228 TNamed(rec),
229
230 fUniformField(rec.fUniformField),
231 fRunVertexFinder(rec.fRunVertexFinder),
232 fRunHLTTracking(rec.fRunHLTTracking),
233 fRunMuonTracking(rec.fRunMuonTracking),
234 fStopOnError(rec.fStopOnError),
235 fWriteAlignmentData(rec.fWriteAlignmentData),
236 fWriteESDfriend(rec.fWriteESDfriend),
237 fWriteAOD(rec.fWriteAOD),
238 fFillTriggerESD(rec.fFillTriggerESD),
239
240 fRunLocalReconstruction(rec.fRunLocalReconstruction),
241 fRunTracking(rec.fRunTracking),
242 fFillESD(rec.fFillESD),
243 fGAliceFileName(rec.fGAliceFileName),
244 fInput(rec.fInput),
245 fEquipIdMap(rec.fEquipIdMap),
246 fFirstEvent(rec.fFirstEvent),
247 fLastEvent(rec.fLastEvent),
248 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
249 fCheckPointLevel(0),
250 fOptions(),
251 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
252 fLoadAlignData(rec.fLoadAlignData),
253 fESDPar(rec.fESDPar),
254
255 fRunLoader(NULL),
256 fRawReader(NULL),
257
258 fVertexer(NULL),
259 fDiamondProfile(NULL),
260
261 fAlignObjArray(rec.fAlignObjArray),
262 fCDBUri(rec.fCDBUri),
263 fSpecCDBUri()
264{
265// copy constructor
266
267 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
268 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
269 }
270 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
271 fReconstructor[iDet] = NULL;
272 fLoader[iDet] = NULL;
273 fTracker[iDet] = NULL;
274 }
275 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
276 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
277 }
278}
279
280//_____________________________________________________________________________
281AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
282{
283// assignment operator
284
285 this->~AliReconstruction();
286 new(this) AliReconstruction(rec);
287 return *this;
288}
289
290//_____________________________________________________________________________
291AliReconstruction::~AliReconstruction()
292{
293// clean up
294
295 CleanUp();
296 fOptions.Delete();
297 fSpecCDBUri.Delete();
298}
299
300//_____________________________________________________________________________
301void AliReconstruction::InitCDBStorage()
302{
303// activate a default CDB storage
304// First check if we have any CDB storage set, because it is used
305// to retrieve the calibration and alignment constants
306
307 AliCDBManager* man = AliCDBManager::Instance();
308 if (man->IsDefaultStorageSet())
309 {
310 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
311 AliWarning("Default CDB storage has been already set !");
312 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
313 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
314 fCDBUri = "";
315 }
316 else {
317 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
318 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
319 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
320 man->SetDefaultStorage(fCDBUri);
321 }
322
323 // Now activate the detector specific CDB storage locations
324 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
325 TObject* obj = fSpecCDBUri[i];
326 if (!obj) continue;
327 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
328 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
329 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
330 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
331 }
332 man->Print();
333}
334
335//_____________________________________________________________________________
336void AliReconstruction::SetDefaultStorage(const char* uri) {
337// Store the desired default CDB storage location
338// Activate it later within the Run() method
339
340 fCDBUri = uri;
341
342}
343
344//_____________________________________________________________________________
345void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
346// Store a detector-specific CDB storage location
347// Activate it later within the Run() method
348
349 AliCDBPath aPath(calibType);
350 if(!aPath.IsValid()){
351 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
352 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
353 if(!strcmp(calibType, fgkDetectorName[iDet])) {
354 aPath.SetPath(Form("%s/*", calibType));
355 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
356 break;
357 }
358 }
359 if(!aPath.IsValid()){
360 AliError(Form("Not a valid path or detector: %s", calibType));
361 return;
362 }
363 }
364
365 // check that calibType refers to a "valid" detector name
366 Bool_t isDetector = kFALSE;
367 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
368 TString detName = fgkDetectorName[iDet];
369 if(aPath.GetLevel0() == detName) {
370 isDetector = kTRUE;
371 break;
372 }
373 }
374
375 if(!isDetector) {
376 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
377 return;
378 }
379
380 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
381 if (obj) fSpecCDBUri.Remove(obj);
382 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
383
384}
385
386
387
388
389//_____________________________________________________________________________
390Bool_t AliReconstruction::SetRunNumber()
391{
392 // The method is called in Run() in order
393 // to set a correct run number.
394 // In case of raw data reconstruction the
395 // run number is taken from the raw data header
396
397 if(AliCDBManager::Instance()->GetRun() < 0) {
398 if (!fRunLoader) {
399 AliError("No run loader is found !");
400 return kFALSE;
401 }
402 // read run number from gAlice
403 if(fRunLoader->GetAliRun())
404 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
405 else {
406 if(fRawReader) {
407 if(fRawReader->NextEvent()) {
408 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
409 fRawReader->RewindEvents();
410 }
411 else {
412 AliError("No raw-data events found !");
413 return kFALSE;
414 }
415 }
416 else {
417 AliError("Neither gAlice nor RawReader objects are found !");
418 return kFALSE;
419 }
420 }
421 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
422 }
423 return kTRUE;
424}
425
426//_____________________________________________________________________________
427Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
428{
429 // Read the alignment objects from CDB.
430 // Each detector is supposed to have the
431 // alignment objects in DET/Align/Data CDB path.
432 // All the detector objects are then collected,
433 // sorted by geometry level (starting from ALIC) and
434 // then applied to the TGeo geometry.
435 // Finally an overlaps check is performed.
436
437 // Load alignment data from CDB and fill fAlignObjArray
438 if(fLoadAlignFromCDB){
439
440 TString detStr = detectors;
441 TString loadAlObjsListOfDets = "";
442
443 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
444 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
445 loadAlObjsListOfDets += fgkDetectorName[iDet];
446 loadAlObjsListOfDets += " ";
447 } // end loop over detectors
448 AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
449 }else{
450 // Check if the array with alignment objects was
451 // provided by the user. If yes, apply the objects
452 // to the present TGeo geometry
453 if (fAlignObjArray) {
454 if (gGeoManager && gGeoManager->IsClosed()) {
455 if (AliGeomManager::ApplyAlignObjsToGeom(*fAlignObjArray) == kFALSE) {
456 AliError("The misalignment of one or more volumes failed!"
457 "Compare the list of simulated detectors and the list of detector alignment data!");
458 return kFALSE;
459 }
460 }
461 else {
462 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
463 return kFALSE;
464 }
465 }
466 }
467
468 delete fAlignObjArray; fAlignObjArray=0;
469
470 return kTRUE;
471}
472
473//_____________________________________________________________________________
474void AliReconstruction::SetGAliceFile(const char* fileName)
475{
476// set the name of the galice file
477
478 fGAliceFileName = fileName;
479}
480
481//_____________________________________________________________________________
482void AliReconstruction::SetOption(const char* detector, const char* option)
483{
484// set options for the reconstruction of a detector
485
486 TObject* obj = fOptions.FindObject(detector);
487 if (obj) fOptions.Remove(obj);
488 fOptions.Add(new TNamed(detector, option));
489}
490
491
492//_____________________________________________________________________________
493Bool_t AliReconstruction::Run(const char* input)
494{
495// run the reconstruction
496
497 // set the input
498 if (!input) input = fInput.Data();
499 TString fileName(input);
500 if (fileName.EndsWith("/")) {
501 fRawReader = new AliRawReaderFile(fileName);
502 } else if (fileName.EndsWith(".root")) {
503 fRawReader = new AliRawReaderRoot(fileName);
504 } else if (!fileName.IsNull()) {
505 fRawReader = new AliRawReaderDate(fileName);
506 fRawReader->SelectEvents(7);
507 }
508 if (!fEquipIdMap.IsNull() && fRawReader)
509 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
510
511
512 // get the run loader
513 if (!InitRunLoader()) return kFALSE;
514
515 // Initialize the CDB storage
516 InitCDBStorage();
517
518 // Set run number in CDBManager (if it is not already set by the user)
519 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
520
521 // Import ideal TGeo geometry and apply misalignment
522 if (!gGeoManager) {
523 TString geom(gSystem->DirName(fGAliceFileName));
524 geom += "/geometry.root";
525 AliGeomManager::LoadGeometry(geom.Data());
526 if (!gGeoManager) if (fStopOnError) return kFALSE;
527 }
528
529 AliCDBManager* man = AliCDBManager::Instance();
530 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
531
532 // local reconstruction
533 if (!fRunLocalReconstruction.IsNull()) {
534 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
535 if (fStopOnError) {CleanUp(); return kFALSE;}
536 }
537 }
538// if (!fRunVertexFinder && fRunTracking.IsNull() &&
539// fFillESD.IsNull()) return kTRUE;
540
541 // get vertexer
542 if (fRunVertexFinder && !CreateVertexer()) {
543 if (fStopOnError) {
544 CleanUp();
545 return kFALSE;
546 }
547 }
548
549 // get trackers
550 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
551 if (fStopOnError) {
552 CleanUp();
553 return kFALSE;
554 }
555 }
556
557
558 TStopwatch stopwatch;
559 stopwatch.Start();
560
561 // get the possibly already existing ESD file and tree
562 AliESD* esd = new AliESD(); AliESD* hltesd = new AliESD();
563 TFile* fileOld = NULL;
564 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
565 if (!gSystem->AccessPathName("AliESDs.root")){
566 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
567 fileOld = TFile::Open("AliESDs.old.root");
568 if (fileOld && fileOld->IsOpen()) {
569 treeOld = (TTree*) fileOld->Get("esdTree");
570 if (treeOld)esd->ReadFromTree(treeOld);
571 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
572 if (hlttreeOld) hltesd->ReadFromTree(hlttreeOld);
573 }
574 }
575
576 // create the ESD output file and tree
577 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
578 file->SetCompressionLevel(2);
579 if (!file->IsOpen()) {
580 AliError("opening AliESDs.root failed");
581 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
582 }
583
584 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
585 esd = new AliESD();
586 esd->CreateStdContent();
587 esd->WriteToTree(tree);
588
589 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
590 hltesd = new AliESD();
591 hltesd->CreateStdContent();
592 hltesd->WriteToTree(hlttree);
593
594 /* CKB Why?
595 delete esd; delete hltesd;
596 esd = NULL; hltesd = NULL;
597 */
598 // create the branch with ESD additions
599 AliESDfriend *esdf = 0;
600 if (fWriteESDfriend) {
601 esdf = new AliESDfriend();
602 TBranch *br=tree->Branch("ESDfriend.","AliESDfriend", &esdf);
603 br->SetFile("AliESDfriends.root");
604 esd->AddObject(esdf);
605 }
606
607 // Get the diamond profile from OCDB
608 AliCDBEntry* entry = AliCDBManager::Instance()
609 ->Get("GRP/Calib/MeanVertex");
610
611 if(entry) {
612 fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
613 } else {
614 AliError("No diamond profile found in OCDB!");
615 }
616
617 AliVertexerTracks tVertexer(AliTracker::GetBz());
618 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
619
620 // loop over events
621 if (fRawReader) fRawReader->RewindEvents();
622
623 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
624 if (fRawReader) fRawReader->NextEvent();
625 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
626 // copy old ESD to the new one
627 if (treeOld) {
628 esd->ReadFromTree(treeOld);
629 treeOld->GetEntry(iEvent);
630 }
631 tree->Fill();
632 if (hlttreeOld) {
633 esd->ReadFromTree(hlttreeOld);
634 hlttreeOld->GetEntry(iEvent);
635 }
636 hlttree->Fill();
637 continue;
638 }
639
640 AliInfo(Form("processing event %d", iEvent));
641 fRunLoader->GetEvent(iEvent);
642
643 char aFileName[256];
644 sprintf(aFileName, "ESD_%d.%d_final.root",
645 fRunLoader->GetHeader()->GetRun(),
646 fRunLoader->GetHeader()->GetEventNrInRun());
647 if (!gSystem->AccessPathName(aFileName)) continue;
648
649 // local reconstruction
650 if (!fRunLocalReconstruction.IsNull()) {
651 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
652 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
653 }
654 }
655
656
657 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
658 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
659 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
660 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
661
662 // Set magnetic field from the tracker
663 esd->SetMagneticField(AliTracker::GetBz());
664 hltesd->SetMagneticField(AliTracker::GetBz());
665
666
667
668 // Fill raw-data error log into the ESD
669 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
670
671 // vertex finder
672 if (fRunVertexFinder) {
673 if (!ReadESD(esd, "vertex")) {
674 if (!RunVertexFinder(esd)) {
675 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
676 }
677 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
678 }
679 }
680
681 // HLT tracking
682 if (!fRunTracking.IsNull()) {
683 if (fRunHLTTracking) {
684 hltesd->SetVertex(esd->GetVertex());
685 if (!RunHLTTracking(hltesd)) {
686 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
687 }
688 }
689 }
690
691 // Muon tracking
692 if (!fRunTracking.IsNull()) {
693 if (fRunMuonTracking) {
694 if (!RunMuonTracking(esd)) {
695 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
696 }
697 }
698 }
699
700 // barrel tracking
701 if (!fRunTracking.IsNull()) {
702 if (!ReadESD(esd, "tracking")) {
703 if (!RunTracking(esd)) {
704 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
705 }
706 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
707 }
708 }
709
710 // fill ESD
711 if (!fFillESD.IsNull()) {
712 if (!FillESD(esd, fFillESD)) {
713 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
714 }
715 }
716 // fill Event header information from the RawEventHeader
717 if (fRawReader){FillRawEventHeaderESD(esd);}
718
719 // combined PID
720 AliESDpid::MakePID(esd);
721 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
722
723 if (fFillTriggerESD) {
724 if (!ReadESD(esd, "trigger")) {
725 if (!FillTriggerESD(esd)) {
726 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
727 }
728 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
729 }
730 }
731
732 //Try to improve the reconstructed primary vertex position using the tracks
733 AliESDVertex *pvtx=0;
734 Bool_t dovertex=kTRUE;
735 TObject* obj = fOptions.FindObject("ITS");
736 if (obj) {
737 TString optITS = obj->GetTitle();
738 if (optITS.Contains("cosmics") || optITS.Contains("COSMICS"))
739 dovertex=kFALSE;
740 }
741 if(dovertex) pvtx=tVertexer.FindPrimaryVertex(esd);
742 if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
743
744 if (pvtx)
745 if (pvtx->GetStatus()) {
746 // Store the improved primary vertex
747 esd->SetPrimaryVertex(pvtx);
748 // Propagate the tracks to the DCA to the improved primary vertex
749 Double_t somethingbig = 777.;
750 Double_t bz = esd->GetMagneticField();
751 Int_t nt=esd->GetNumberOfTracks();
752 while (nt--) {
753 AliESDtrack *t = esd->GetTrack(nt);
754 t->RelateToVertex(pvtx, bz, somethingbig);
755 }
756 }
757
758 {
759 // V0 finding
760 AliV0vertexer vtxer;
761 vtxer.Tracks2V0vertices(esd);
762
763 // Cascade finding
764 AliCascadeVertexer cvtxer;
765 cvtxer.V0sTracks2CascadeVertices(esd);
766 }
767
768 // write ESD
769 if (fWriteESDfriend) {
770 new (esdf) AliESDfriend(); // Reset...
771 esd->GetESDfriend(esdf);
772 }
773 tree->Fill();
774
775 // write HLT ESD
776 hlttree->Fill();
777
778 if (fCheckPointLevel > 0) WriteESD(esd, "final");
779 esd->Reset();
780 hltesd->Reset();
781 // esdf->Reset();
782 // delete esdf; esdf = 0;
783 }
784
785
786 tree->GetUserInfo()->Add(esd);
787 hlttree->GetUserInfo()->Add(hltesd);
788
789
790
791 if(fESDPar.Contains("ESD.par")){
792 AliInfo("Attaching ESD.par to Tree");
793 TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
794 tree->GetUserInfo()->Add(fn);
795 }
796
797
798 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
799 stopwatch.RealTime(),stopwatch.CpuTime()));
800
801 file->cd();
802 if (fWriteESDfriend)
803 tree->SetBranchStatus("ESDfriend*",0);
804 tree->Write();
805 hlttree->Write();
806
807 if (fWriteAOD) {
808 TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
809 ESDFile2AODFile(file, aodFile);
810 aodFile->Close();
811 }
812
813 gROOT->cd();
814 // Create tags for the events in the ESD tree (the ESD tree is always present)
815 // In case of empty events the tags will contain dummy values
816 CreateTag(file);
817
818 CleanUp(file, fileOld);
819
820 return kTRUE;
821}
822
823
824//_____________________________________________________________________________
825Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
826{
827// run the local reconstruction
828
829 TStopwatch stopwatch;
830 stopwatch.Start();
831
832 AliCDBManager* man = AliCDBManager::Instance();
833 Bool_t origCache = man->GetCacheFlag();
834
835 TString detStr = detectors;
836 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
837 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
838 AliReconstructor* reconstructor = GetReconstructor(iDet);
839 if (!reconstructor) continue;
840 if (reconstructor->HasLocalReconstruction()) continue;
841
842 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
843 TStopwatch stopwatchDet;
844 stopwatchDet.Start();
845
846 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
847
848 man->SetCacheFlag(kTRUE);
849 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
850 man->GetAll(calibPath); // entries are cached!
851
852 if (fRawReader) {
853 fRawReader->RewindEvents();
854 reconstructor->Reconstruct(fRunLoader, fRawReader);
855 } else {
856 reconstructor->Reconstruct(fRunLoader);
857 }
858 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
859 fgkDetectorName[iDet],
860 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
861
862 // unload calibration data
863 man->UnloadFromCache(calibPath);
864 //man->ClearCache();
865 }
866
867 man->SetCacheFlag(origCache);
868
869 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
870 AliError(Form("the following detectors were not found: %s",
871 detStr.Data()));
872 if (fStopOnError) return kFALSE;
873 }
874
875 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
876 stopwatch.RealTime(),stopwatch.CpuTime()));
877
878 return kTRUE;
879}
880
881//_____________________________________________________________________________
882Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
883{
884// run the local reconstruction
885
886 TStopwatch stopwatch;
887 stopwatch.Start();
888
889 TString detStr = detectors;
890 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
891 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
892 AliReconstructor* reconstructor = GetReconstructor(iDet);
893 if (!reconstructor) continue;
894 AliLoader* loader = fLoader[iDet];
895
896 // conversion of digits
897 if (fRawReader && reconstructor->HasDigitConversion()) {
898 AliInfo(Form("converting raw data digits into root objects for %s",
899 fgkDetectorName[iDet]));
900 TStopwatch stopwatchDet;
901 stopwatchDet.Start();
902 loader->LoadDigits("update");
903 loader->CleanDigits();
904 loader->MakeDigitsContainer();
905 TTree* digitsTree = loader->TreeD();
906 reconstructor->ConvertDigits(fRawReader, digitsTree);
907 loader->WriteDigits("OVERWRITE");
908 loader->UnloadDigits();
909 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
910 fgkDetectorName[iDet],
911 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
912 }
913
914 // local reconstruction
915 if (!reconstructor->HasLocalReconstruction()) continue;
916 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
917 TStopwatch stopwatchDet;
918 stopwatchDet.Start();
919 loader->LoadRecPoints("update");
920 loader->CleanRecPoints();
921 loader->MakeRecPointsContainer();
922 TTree* clustersTree = loader->TreeR();
923 if (fRawReader && !reconstructor->HasDigitConversion()) {
924 reconstructor->Reconstruct(fRawReader, clustersTree);
925 } else {
926 loader->LoadDigits("read");
927 TTree* digitsTree = loader->TreeD();
928 if (!digitsTree) {
929 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
930 if (fStopOnError) return kFALSE;
931 } else {
932 reconstructor->Reconstruct(digitsTree, clustersTree);
933 }
934 loader->UnloadDigits();
935 }
936 loader->WriteRecPoints("OVERWRITE");
937 loader->UnloadRecPoints();
938 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
939 fgkDetectorName[iDet],
940 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
941 }
942
943 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
944 AliError(Form("the following detectors were not found: %s",
945 detStr.Data()));
946 if (fStopOnError) return kFALSE;
947 }
948
949 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
950 stopwatch.RealTime(),stopwatch.CpuTime()));
951
952 return kTRUE;
953}
954
955//_____________________________________________________________________________
956Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
957{
958// run the barrel tracking
959
960 TStopwatch stopwatch;
961 stopwatch.Start();
962
963 AliESDVertex* vertex = NULL;
964 Double_t vtxPos[3] = {0, 0, 0};
965 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
966 TArrayF mcVertex(3);
967 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
968 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
969 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
970 }
971
972 if (fVertexer) {
973 if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
974 AliInfo("running the ITS vertex finder");
975 if (fLoader[0]) fLoader[0]->LoadRecPoints();
976 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
977 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
978 if(!vertex){
979 AliWarning("Vertex not found");
980 vertex = new AliESDVertex();
981 vertex->SetName("default");
982 }
983 else {
984 vertex->SetName("reconstructed");
985 }
986
987 } else {
988 AliInfo("getting the primary vertex from MC");
989 vertex = new AliESDVertex(vtxPos, vtxErr);
990 }
991
992 if (vertex) {
993 vertex->GetXYZ(vtxPos);
994 vertex->GetSigmaXYZ(vtxErr);
995 } else {
996 AliWarning("no vertex reconstructed");
997 vertex = new AliESDVertex(vtxPos, vtxErr);
998 }
999 esd->SetVertex(vertex);
1000 // if SPD multiplicity has been determined, it is stored in the ESD
1001 AliMultiplicity *mult = fVertexer->GetMultiplicity();
1002 if(mult)esd->SetMultiplicity(mult);
1003
1004 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1005 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1006 }
1007 delete vertex;
1008
1009 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1010 stopwatch.RealTime(),stopwatch.CpuTime()));
1011
1012 return kTRUE;
1013}
1014
1015//_____________________________________________________________________________
1016Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1017{
1018// run the HLT barrel tracking
1019
1020 TStopwatch stopwatch;
1021 stopwatch.Start();
1022
1023 if (!fRunLoader) {
1024 AliError("Missing runLoader!");
1025 return kFALSE;
1026 }
1027
1028 AliInfo("running HLT tracking");
1029
1030 // Get a pointer to the HLT reconstructor
1031 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1032 if (!reconstructor) return kFALSE;
1033
1034 // TPC + ITS
1035 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1036 TString detName = fgkDetectorName[iDet];
1037 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1038 reconstructor->SetOption(detName.Data());
1039 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1040 if (!tracker) {
1041 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1042 if (fStopOnError) return kFALSE;
1043 continue;
1044 }
1045 Double_t vtxPos[3];
1046 Double_t vtxErr[3]={0.005,0.005,0.010};
1047 const AliESDVertex *vertex = esd->GetVertex();
1048 vertex->GetXYZ(vtxPos);
1049 tracker->SetVertex(vtxPos,vtxErr);
1050 if(iDet != 1) {
1051 fLoader[iDet]->LoadRecPoints("read");
1052 TTree* tree = fLoader[iDet]->TreeR();
1053 if (!tree) {
1054 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1055 return kFALSE;
1056 }
1057 tracker->LoadClusters(tree);
1058 }
1059 if (tracker->Clusters2Tracks(esd) != 0) {
1060 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1061 return kFALSE;
1062 }
1063 if(iDet != 1) {
1064 tracker->UnloadClusters();
1065 }
1066 delete tracker;
1067 }
1068
1069 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1070 stopwatch.RealTime(),stopwatch.CpuTime()));
1071
1072 return kTRUE;
1073}
1074
1075//_____________________________________________________________________________
1076Bool_t AliReconstruction::RunMuonTracking(AliESD*& esd)
1077{
1078// run the muon spectrometer tracking
1079
1080 TStopwatch stopwatch;
1081 stopwatch.Start();
1082
1083 if (!fRunLoader) {
1084 AliError("Missing runLoader!");
1085 return kFALSE;
1086 }
1087 Int_t iDet = 7; // for MUON
1088
1089 AliInfo("is running...");
1090
1091 // Get a pointer to the MUON reconstructor
1092 AliReconstructor *reconstructor = GetReconstructor(iDet);
1093 if (!reconstructor) return kFALSE;
1094
1095
1096 TString detName = fgkDetectorName[iDet];
1097 AliDebug(1, Form("%s tracking", detName.Data()));
1098 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1099 if (!tracker) {
1100 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1101 return kFALSE;
1102 }
1103
1104 // create Tracks
1105 fLoader[iDet]->LoadTracks("update");
1106 fLoader[iDet]->CleanTracks();
1107 fLoader[iDet]->MakeTracksContainer();
1108
1109 // read RecPoints
1110 fLoader[iDet]->LoadRecPoints("read");
1111 tracker->LoadClusters(fLoader[iDet]->TreeR());
1112
1113 Int_t rv = tracker->Clusters2Tracks(esd);
1114
1115 fLoader[iDet]->UnloadRecPoints();
1116
1117 if ( rv )
1118 {
1119 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1120 return kFALSE;
1121 }
1122
1123 tracker->UnloadClusters();
1124
1125 fLoader[iDet]->UnloadRecPoints();
1126
1127 fLoader[iDet]->WriteTracks("OVERWRITE");
1128 fLoader[iDet]->UnloadTracks();
1129
1130 delete tracker;
1131
1132 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1133 stopwatch.RealTime(),stopwatch.CpuTime()));
1134
1135 return kTRUE;
1136}
1137
1138
1139//_____________________________________________________________________________
1140Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1141{
1142// run the barrel tracking
1143
1144 TStopwatch stopwatch;
1145 stopwatch.Start();
1146
1147 AliInfo("running tracking");
1148
1149 //Fill the ESD with the T0 info (will be used by the TOF)
1150 if (fReconstructor[11])
1151 GetReconstructor(11)->FillESD(fRunLoader, esd);
1152
1153 // pass 1: TPC + ITS inwards
1154 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1155 if (!fTracker[iDet]) continue;
1156 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1157
1158 // load clusters
1159 fLoader[iDet]->LoadRecPoints("read");
1160 TTree* tree = fLoader[iDet]->TreeR();
1161 if (!tree) {
1162 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1163 return kFALSE;
1164 }
1165 fTracker[iDet]->LoadClusters(tree);
1166
1167 // run tracking
1168 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1169 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1170 return kFALSE;
1171 }
1172 if (fCheckPointLevel > 1) {
1173 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1174 }
1175 // preliminary PID in TPC needed by the ITS tracker
1176 if (iDet == 1) {
1177 GetReconstructor(1)->FillESD(fRunLoader, esd);
1178 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1179 AliESDpid::MakePID(esd);
1180 }
1181 }
1182
1183 // pass 2: ALL backwards
1184 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1185 if (!fTracker[iDet]) continue;
1186 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1187
1188 // load clusters
1189 if (iDet > 1) { // all except ITS, TPC
1190 TTree* tree = NULL;
1191 fLoader[iDet]->LoadRecPoints("read");
1192 tree = fLoader[iDet]->TreeR();
1193 if (!tree) {
1194 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1195 return kFALSE;
1196 }
1197 fTracker[iDet]->LoadClusters(tree);
1198 }
1199
1200 // run tracking
1201 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1202 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1203 // return kFALSE;
1204 }
1205 if (fCheckPointLevel > 1) {
1206 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1207 }
1208
1209 // unload clusters
1210 if (iDet > 2) { // all except ITS, TPC, TRD
1211 fTracker[iDet]->UnloadClusters();
1212 fLoader[iDet]->UnloadRecPoints();
1213 }
1214 // updated PID in TPC needed by the ITS tracker -MI
1215 if (iDet == 1) {
1216 GetReconstructor(1)->FillESD(fRunLoader, esd);
1217 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1218 AliESDpid::MakePID(esd);
1219 }
1220 }
1221
1222 // write space-points to the ESD in case alignment data output
1223 // is switched on
1224 if (fWriteAlignmentData)
1225 WriteAlignmentData(esd);
1226
1227 // pass 3: TRD + TPC + ITS refit inwards
1228 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1229 if (!fTracker[iDet]) continue;
1230 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1231
1232 // run tracking
1233 if (fTracker[iDet]->RefitInward(esd) != 0) {
1234 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1235 // return kFALSE;
1236 }
1237 if (fCheckPointLevel > 1) {
1238 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1239 }
1240
1241 // unload clusters
1242 fTracker[iDet]->UnloadClusters();
1243 fLoader[iDet]->UnloadRecPoints();
1244 }
1245 //
1246 // Propagate track to the vertex - if not done by ITS
1247 //
1248 Int_t ntracks = esd->GetNumberOfTracks();
1249 for (Int_t itrack=0; itrack<ntracks; itrack++){
1250 const Double_t kRadius = 3; // beam pipe radius
1251 const Double_t kMaxStep = 5; // max step
1252 const Double_t kMaxD = 123456; // max distance to prim vertex
1253 Double_t fieldZ = AliTracker::GetBz(); //
1254 AliESDtrack * track = esd->GetTrack(itrack);
1255 if (!track) continue;
1256 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1257 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1258 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1259 }
1260
1261 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1262 stopwatch.RealTime(),stopwatch.CpuTime()));
1263
1264 return kTRUE;
1265}
1266
1267//_____________________________________________________________________________
1268Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1269{
1270// fill the event summary data
1271
1272 TStopwatch stopwatch;
1273 stopwatch.Start();
1274 AliInfo("filling ESD");
1275
1276 TString detStr = detectors;
1277 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1278 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1279 AliReconstructor* reconstructor = GetReconstructor(iDet);
1280 if (!reconstructor) continue;
1281
1282 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1283 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1284 TTree* clustersTree = NULL;
1285 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1286 fLoader[iDet]->LoadRecPoints("read");
1287 clustersTree = fLoader[iDet]->TreeR();
1288 if (!clustersTree) {
1289 AliError(Form("Can't get the %s clusters tree",
1290 fgkDetectorName[iDet]));
1291 if (fStopOnError) return kFALSE;
1292 }
1293 }
1294 if (fRawReader && !reconstructor->HasDigitConversion()) {
1295 reconstructor->FillESD(fRawReader, clustersTree, esd);
1296 } else {
1297 TTree* digitsTree = NULL;
1298 if (fLoader[iDet]) {
1299 fLoader[iDet]->LoadDigits("read");
1300 digitsTree = fLoader[iDet]->TreeD();
1301 if (!digitsTree) {
1302 AliError(Form("Can't get the %s digits tree",
1303 fgkDetectorName[iDet]));
1304 if (fStopOnError) return kFALSE;
1305 }
1306 }
1307 reconstructor->FillESD(digitsTree, clustersTree, esd);
1308 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1309 }
1310 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1311 fLoader[iDet]->UnloadRecPoints();
1312 }
1313
1314 if (fRawReader) {
1315 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1316 } else {
1317 reconstructor->FillESD(fRunLoader, esd);
1318 }
1319 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1320 }
1321 }
1322
1323 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1324 AliError(Form("the following detectors were not found: %s",
1325 detStr.Data()));
1326 if (fStopOnError) return kFALSE;
1327 }
1328
1329 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1330 stopwatch.RealTime(),stopwatch.CpuTime()));
1331
1332 return kTRUE;
1333}
1334
1335//_____________________________________________________________________________
1336Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1337{
1338 // Reads the trigger decision which is
1339 // stored in Trigger.root file and fills
1340 // the corresponding esd entries
1341
1342 AliInfo("Filling trigger information into the ESD");
1343
1344 if (fRawReader) {
1345 AliCTPRawStream input(fRawReader);
1346 if (!input.Next()) {
1347 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1348 return kFALSE;
1349 }
1350 esd->SetTriggerMask(input.GetClassMask());
1351 esd->SetTriggerCluster(input.GetClusterMask());
1352 }
1353 else {
1354 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1355 if (runloader) {
1356 if (!runloader->LoadTrigger()) {
1357 AliCentralTrigger *aCTP = runloader->GetTrigger();
1358 esd->SetTriggerMask(aCTP->GetClassMask());
1359 esd->SetTriggerCluster(aCTP->GetClusterMask());
1360 }
1361 else {
1362 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1363 return kFALSE;
1364 }
1365 }
1366 else {
1367 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1368 return kFALSE;
1369 }
1370 }
1371
1372 return kTRUE;
1373}
1374
1375
1376
1377
1378
1379//_____________________________________________________________________________
1380Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1381{
1382 //
1383 // Filling information from RawReader Header
1384 //
1385
1386 AliInfo("Filling information from RawReader Header");
1387 esd->SetBunchCrossNumber(0);
1388 esd->SetOrbitNumber(0);
1389 esd->SetPeriodNumber(0);
1390 esd->SetTimeStamp(0);
1391 esd->SetEventType(0);
1392 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1393 if (eventHeader){
1394
1395 const UInt_t *id = eventHeader->GetP("Id");
1396 esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1397 esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1398 esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1399
1400 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1401 esd->SetEventType((eventHeader->Get("Type")));
1402 }
1403
1404 return kTRUE;
1405}
1406
1407
1408//_____________________________________________________________________________
1409Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1410{
1411// check whether detName is contained in detectors
1412// if yes, it is removed from detectors
1413
1414 // check if all detectors are selected
1415 if ((detectors.CompareTo("ALL") == 0) ||
1416 detectors.BeginsWith("ALL ") ||
1417 detectors.EndsWith(" ALL") ||
1418 detectors.Contains(" ALL ")) {
1419 detectors = "ALL";
1420 return kTRUE;
1421 }
1422
1423 // search for the given detector
1424 Bool_t result = kFALSE;
1425 if ((detectors.CompareTo(detName) == 0) ||
1426 detectors.BeginsWith(detName+" ") ||
1427 detectors.EndsWith(" "+detName) ||
1428 detectors.Contains(" "+detName+" ")) {
1429 detectors.ReplaceAll(detName, "");
1430 result = kTRUE;
1431 }
1432
1433 // clean up the detectors string
1434 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1435 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1436 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1437
1438 return result;
1439}
1440
1441//_____________________________________________________________________________
1442Bool_t AliReconstruction::InitRunLoader()
1443{
1444// get or create the run loader
1445
1446 if (gAlice) delete gAlice;
1447 gAlice = NULL;
1448
1449 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1450 // load all base libraries to get the loader classes
1451 TString libs = gSystem->GetLibraries();
1452 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1453 TString detName = fgkDetectorName[iDet];
1454 if (detName == "HLT") continue;
1455 if (libs.Contains("lib" + detName + "base.so")) continue;
1456 gSystem->Load("lib" + detName + "base.so");
1457 }
1458 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1459 if (!fRunLoader) {
1460 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1461 CleanUp();
1462 return kFALSE;
1463 }
1464 fRunLoader->CdGAFile();
1465 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1466 if (fRunLoader->LoadgAlice() == 0) {
1467 gAlice = fRunLoader->GetAliRun();
1468 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1469 }
1470 }
1471 if (!gAlice && !fRawReader) {
1472 AliError(Form("no gAlice object found in file %s",
1473 fGAliceFileName.Data()));
1474 CleanUp();
1475 return kFALSE;
1476 }
1477
1478 } else { // galice.root does not exist
1479 if (!fRawReader) {
1480 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1481 CleanUp();
1482 return kFALSE;
1483 }
1484 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1485 AliConfig::GetDefaultEventFolderName(),
1486 "recreate");
1487 if (!fRunLoader) {
1488 AliError(Form("could not create run loader in file %s",
1489 fGAliceFileName.Data()));
1490 CleanUp();
1491 return kFALSE;
1492 }
1493 fRunLoader->MakeTree("E");
1494 Int_t iEvent = 0;
1495 while (fRawReader->NextEvent()) {
1496 fRunLoader->SetEventNumber(iEvent);
1497 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1498 iEvent, iEvent);
1499 fRunLoader->MakeTree("H");
1500 fRunLoader->TreeE()->Fill();
1501 iEvent++;
1502 }
1503 fRawReader->RewindEvents();
1504 if (fNumberOfEventsPerFile > 0)
1505 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1506 else
1507 fRunLoader->SetNumberOfEventsPerFile(iEvent);
1508 fRunLoader->WriteHeader("OVERWRITE");
1509 fRunLoader->CdGAFile();
1510 fRunLoader->Write(0, TObject::kOverwrite);
1511// AliTracker::SetFieldMap(???);
1512 }
1513
1514 return kTRUE;
1515}
1516
1517//_____________________________________________________________________________
1518AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1519{
1520// get the reconstructor object and the loader for a detector
1521
1522 if (fReconstructor[iDet]) return fReconstructor[iDet];
1523
1524 // load the reconstructor object
1525 TPluginManager* pluginManager = gROOT->GetPluginManager();
1526 TString detName = fgkDetectorName[iDet];
1527 TString recName = "Ali" + detName + "Reconstructor";
1528 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1529
1530 if (detName == "HLT") {
1531 if (!gROOT->GetClass("AliLevel3")) {
1532 gSystem->Load("libAliHLTSrc.so");
1533 gSystem->Load("libAliHLTMisc.so");
1534 gSystem->Load("libAliHLTHough.so");
1535 gSystem->Load("libAliHLTComp.so");
1536 }
1537 }
1538
1539 AliReconstructor* reconstructor = NULL;
1540 // first check if a plugin is defined for the reconstructor
1541 TPluginHandler* pluginHandler =
1542 pluginManager->FindHandler("AliReconstructor", detName);
1543 // if not, add a plugin for it
1544 if (!pluginHandler) {
1545 AliDebug(1, Form("defining plugin for %s", recName.Data()));
1546 TString libs = gSystem->GetLibraries();
1547 if (libs.Contains("lib" + detName + "base.so") ||
1548 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1549 pluginManager->AddHandler("AliReconstructor", detName,
1550 recName, detName + "rec", recName + "()");
1551 } else {
1552 pluginManager->AddHandler("AliReconstructor", detName,
1553 recName, detName, recName + "()");
1554 }
1555 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1556 }
1557 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1558 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1559 }
1560 if (reconstructor) {
1561 TObject* obj = fOptions.FindObject(detName.Data());
1562 if (obj) reconstructor->SetOption(obj->GetTitle());
1563 reconstructor->Init(fRunLoader);
1564 fReconstructor[iDet] = reconstructor;
1565 }
1566
1567 // get or create the loader
1568 if (detName != "HLT") {
1569 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1570 if (!fLoader[iDet]) {
1571 AliConfig::Instance()
1572 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1573 detName, detName);
1574 // first check if a plugin is defined for the loader
1575 pluginHandler =
1576 pluginManager->FindHandler("AliLoader", detName);
1577 // if not, add a plugin for it
1578 if (!pluginHandler) {
1579 TString loaderName = "Ali" + detName + "Loader";
1580 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1581 pluginManager->AddHandler("AliLoader", detName,
1582 loaderName, detName + "base",
1583 loaderName + "(const char*, TFolder*)");
1584 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1585 }
1586 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1587 fLoader[iDet] =
1588 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1589 fRunLoader->GetEventFolder());
1590 }
1591 if (!fLoader[iDet]) { // use default loader
1592 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1593 }
1594 if (!fLoader[iDet]) {
1595 AliWarning(Form("couldn't get loader for %s", detName.Data()));
1596 if (fStopOnError) return NULL;
1597 } else {
1598 fRunLoader->AddLoader(fLoader[iDet]);
1599 fRunLoader->CdGAFile();
1600 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1601 fRunLoader->Write(0, TObject::kOverwrite);
1602 }
1603 }
1604 }
1605
1606 return reconstructor;
1607}
1608
1609//_____________________________________________________________________________
1610Bool_t AliReconstruction::CreateVertexer()
1611{
1612// create the vertexer
1613
1614 fVertexer = NULL;
1615 AliReconstructor* itsReconstructor = GetReconstructor(0);
1616 if (itsReconstructor) {
1617 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1618 }
1619 if (!fVertexer) {
1620 AliWarning("couldn't create a vertexer for ITS");
1621 if (fStopOnError) return kFALSE;
1622 }
1623
1624 return kTRUE;
1625}
1626
1627//_____________________________________________________________________________
1628Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1629{
1630// create the trackers
1631
1632 TString detStr = detectors;
1633 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1634 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1635 AliReconstructor* reconstructor = GetReconstructor(iDet);
1636 if (!reconstructor) continue;
1637 TString detName = fgkDetectorName[iDet];
1638 if (detName == "HLT") {
1639 fRunHLTTracking = kTRUE;
1640 continue;
1641 }
1642 if (detName == "MUON") {
1643 fRunMuonTracking = kTRUE;
1644 continue;
1645 }
1646
1647
1648 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1649 if (!fTracker[iDet] && (iDet < 7)) {
1650 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1651 if (fStopOnError) return kFALSE;
1652 }
1653 }
1654
1655 return kTRUE;
1656}
1657
1658//_____________________________________________________________________________
1659void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1660{
1661// delete trackers and the run loader and close and delete the file
1662
1663 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1664 delete fReconstructor[iDet];
1665 fReconstructor[iDet] = NULL;
1666 fLoader[iDet] = NULL;
1667 delete fTracker[iDet];
1668 fTracker[iDet] = NULL;
1669 }
1670 delete fVertexer;
1671 fVertexer = NULL;
1672 delete fDiamondProfile;
1673 fDiamondProfile = NULL;
1674
1675 delete fRunLoader;
1676 fRunLoader = NULL;
1677 delete fRawReader;
1678 fRawReader = NULL;
1679
1680 if (file) {
1681 file->Close();
1682 delete file;
1683 }
1684
1685 if (fileOld) {
1686 fileOld->Close();
1687 delete fileOld;
1688 gSystem->Unlink("AliESDs.old.root");
1689 }
1690}
1691
1692
1693//_____________________________________________________________________________
1694Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1695{
1696// read the ESD event from a file
1697
1698 if (!esd) return kFALSE;
1699 char fileName[256];
1700 sprintf(fileName, "ESD_%d.%d_%s.root",
1701 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1702 if (gSystem->AccessPathName(fileName)) return kFALSE;
1703
1704 AliInfo(Form("reading ESD from file %s", fileName));
1705 AliDebug(1, Form("reading ESD from file %s", fileName));
1706 TFile* file = TFile::Open(fileName);
1707 if (!file || !file->IsOpen()) {
1708 AliError(Form("opening %s failed", fileName));
1709 delete file;
1710 return kFALSE;
1711 }
1712
1713 gROOT->cd();
1714 delete esd;
1715 esd = (AliESD*) file->Get("ESD");
1716 file->Close();
1717 delete file;
1718 return kTRUE;
1719}
1720
1721//_____________________________________________________________________________
1722void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1723{
1724// write the ESD event to a file
1725
1726 if (!esd) return;
1727 char fileName[256];
1728 sprintf(fileName, "ESD_%d.%d_%s.root",
1729 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1730
1731 AliDebug(1, Form("writing ESD to file %s", fileName));
1732 TFile* file = TFile::Open(fileName, "recreate");
1733 if (!file || !file->IsOpen()) {
1734 AliError(Form("opening %s failed", fileName));
1735 } else {
1736 esd->Write("ESD");
1737 file->Close();
1738 }
1739 delete file;
1740}
1741
1742
1743
1744
1745//_____________________________________________________________________________
1746void AliReconstruction::CreateTag(TFile* file)
1747{
1748 //GRP
1749 Float_t lhcLuminosity = 0.0;
1750 TString lhcState = "test";
1751 UInt_t detectorMask = 0;
1752
1753 /////////////
1754 //muon code//
1755 ////////////
1756 Double_t fMUONMASS = 0.105658369;
1757 //Variables
1758 Double_t fX,fY,fZ ;
1759 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1760 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1761 Int_t fCharge;
1762 TLorentzVector fEPvector;
1763
1764 Float_t fZVertexCut = 10.0;
1765 Float_t fRhoVertexCut = 2.0;
1766
1767 Float_t fLowPtCut = 1.0;
1768 Float_t fHighPtCut = 3.0;
1769 Float_t fVeryHighPtCut = 10.0;
1770 ////////////
1771
1772 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1773
1774 // Creates the tags for all the events in a given ESD file
1775 Int_t ntrack;
1776 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1777 Int_t nPos, nNeg, nNeutr;
1778 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1779 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1780 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1781 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1782 Float_t maxPt = .0, meanPt = .0, totalP = .0;
1783 Int_t fVertexflag;
1784 Int_t iRunNumber = 0;
1785 TString fVertexName("default");
1786
1787 AliRunTag *tag = new AliRunTag();
1788 AliEventTag *evTag = new AliEventTag();
1789 TTree ttag("T","A Tree with event tags");
1790 TBranch * btag = ttag.Branch("AliTAG", &tag);
1791 btag->SetCompressionLevel(9);
1792
1793 AliInfo(Form("Creating the tags......."));
1794
1795 if (!file || !file->IsOpen()) {
1796 AliError(Form("opening failed"));
1797 delete file;
1798 return ;
1799 }
1800 Int_t lastEvent = 0;
1801 TTree *b = (TTree*) file->Get("esdTree");
1802 AliESD *esd = new AliESD();
1803 esd->ReadFromTree(b);
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->SetLHCTag(lhcLuminosity,lhcState);
2012 tag->SetDetectorTag(detectorMask);
2013
2014 tag->SetRunId(iInitRunNumber);
2015 tag->AddEventTag(*evTag);
2016 }
2017 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2018 else lastEvent = fLastEvent;
2019
2020 ttag.Fill();
2021 tag->Clear();
2022
2023 char fileName[256];
2024 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2025 tag->GetRunId(),fFirstEvent,lastEvent );
2026 AliInfo(Form("writing tags to file %s", fileName));
2027 AliDebug(1, Form("writing tags to file %s", fileName));
2028
2029 TFile* ftag = TFile::Open(fileName, "recreate");
2030 ftag->cd();
2031 ttag.Write();
2032 ftag->Close();
2033 file->cd();
2034 delete tag;
2035 delete evTag;
2036}
2037
2038//_____________________________________________________________________________
2039void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
2040{
2041 // write all files from the given esd file to an aod file
2042
2043 // create an AliAOD object
2044 AliAODEvent *aod = new AliAODEvent();
2045 aod->CreateStdContent();
2046
2047 // go to the file
2048 aodFile->cd();
2049
2050 // create the tree
2051 TTree *aodTree = new TTree("AOD", "AliAOD tree");
2052 aodTree->Branch(aod->GetList());
2053
2054 // connect to ESD
2055 TTree *t = (TTree*) esdFile->Get("esdTree");
2056 AliESD *esd = new AliESD();
2057 esd->ReadFromTree(t);
2058
2059 Int_t nEvents = t->GetEntries();
2060
2061 // set arrays and pointers
2062 Float_t posF[3];
2063 Double_t pos[3];
2064 Double_t p[3];
2065 Double_t covVtx[6];
2066 Double_t covTr[21];
2067 Double_t pid[10];
2068
2069 // loop over events and fill them
2070 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2071 t->GetEntry(iEvent);
2072
2073 // Multiplicity information needed by the header (to be revised!)
2074 Int_t nTracks = esd->GetNumberOfTracks();
2075 Int_t nPosTracks = 0;
2076 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2077 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2078
2079 // Update the header
2080 AliAODHeader* header = aod->GetHeader();
2081
2082 header->SetRunNumber (esd->GetRunNumber() );
2083 header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
2084 header->SetOrbitNumber (esd->GetOrbitNumber() );
2085 header->SetPeriodNumber (esd->GetPeriodNumber() );
2086 header->SetTriggerMask (esd->GetTriggerMask() );
2087 header->SetTriggerCluster (esd->GetTriggerCluster() );
2088 header->SetEventType (esd->GetEventType() );
2089 header->SetMagneticField (esd->GetMagneticField() );
2090 header->SetZDCN1Energy (esd->GetZDCN1Energy() );
2091 header->SetZDCP1Energy (esd->GetZDCP1Energy() );
2092 header->SetZDCN2Energy (esd->GetZDCN2Energy() );
2093 header->SetZDCP2Energy (esd->GetZDCP2Energy() );
2094 header->SetZDCEMEnergy (esd->GetZDCEMEnergy() );
2095 header->SetRefMultiplicity (nTracks);
2096 header->SetRefMultiplicityPos(nPosTracks);
2097 header->SetRefMultiplicityNeg(nTracks - nPosTracks);
2098 header->SetMuonMagFieldScale(-999.); // FIXME
2099 header->SetCentrality(-999.); // FIXME
2100//
2101//
2102
2103 Int_t nV0s = esd->GetNumberOfV0s();
2104 Int_t nCascades = esd->GetNumberOfCascades();
2105 Int_t nKinks = esd->GetNumberOfKinks();
2106 Int_t nVertices = nV0s + nCascades + nKinks;
2107
2108 aod->ResetStd(nTracks, nVertices);
2109 AliAODTrack *aodTrack;
2110
2111
2112 // Array to take into account the tracks already added to the AOD
2113 Bool_t * usedTrack = NULL;
2114 if (nTracks>0) {
2115 usedTrack = new Bool_t[nTracks];
2116 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2117 }
2118 // Array to take into account the V0s already added to the AOD
2119 Bool_t * usedV0 = NULL;
2120 if (nV0s>0) {
2121 usedV0 = new Bool_t[nV0s];
2122 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2123 }
2124 // Array to take into account the kinks already added to the AOD
2125 Bool_t * usedKink = NULL;
2126 if (nKinks>0) {
2127 usedKink = new Bool_t[nKinks];
2128 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2129 }
2130
2131 // Access to the AOD container of vertices
2132 TClonesArray &vertices = *(aod->GetVertices());
2133 Int_t jVertices=0;
2134
2135 // Access to the AOD container of tracks
2136 TClonesArray &tracks = *(aod->GetTracks());
2137 Int_t jTracks=0;
2138
2139 // Add primary vertex. The primary tracks will be defined
2140 // after the loops on the composite objects (V0, cascades, kinks)
2141 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2142
2143 vtx->GetXYZ(pos); // position
2144 vtx->GetCovMatrix(covVtx); //covariance matrix
2145
2146 AliAODVertex * primary = new(vertices[jVertices++])
2147 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2148
2149 // Create vertices starting from the most complex objects
2150
2151 // Cascades
2152 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2153 AliESDcascade *cascade = esd->GetCascade(nCascade);
2154
2155 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2156 cascade->GetPosCovXi(covVtx);
2157
2158 // Add the cascade vertex
2159 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2160 covVtx,
2161 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2162 primary,
2163 AliAODVertex::kCascade);
2164
2165 primary->AddDaughter(vcascade);
2166
2167 // Add the V0 from the cascade. The ESD class have to be optimized...
2168 // Now we have to search for the corresponding Vo in the list of V0s
2169 // using the indeces of the positive and negative tracks
2170
2171 Int_t posFromV0 = cascade->GetPindex();
2172 Int_t negFromV0 = cascade->GetNindex();
2173
2174
2175 AliESDv0 * v0 = 0x0;
2176 Int_t indV0 = -1;
2177
2178 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2179
2180 v0 = esd->GetV0(iV0);
2181 Int_t posV0 = v0->GetPindex();
2182 Int_t negV0 = v0->GetNindex();
2183
2184 if (posV0==posFromV0 && negV0==negFromV0) {
2185 indV0 = iV0;
2186 break;
2187 }
2188 }
2189
2190 AliAODVertex * vV0FromCascade = 0x0;
2191
2192 if (indV0>-1 && !usedV0[indV0] ) {
2193
2194 // the V0 exists in the array of V0s and is not used
2195
2196 usedV0[indV0] = kTRUE;
2197
2198 v0->GetXYZ(pos[0], pos[1], pos[2]);
2199 v0->GetPosCov(covVtx);
2200
2201 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2202 covVtx,
2203 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2204 vcascade,
2205 AliAODVertex::kV0);
2206 } else {
2207
2208 // the V0 doesn't exist in the array of V0s or was used
2209 cerr << "Error: event " << iEvent << " cascade " << nCascade
2210 << " The V0 " << indV0
2211 << " doesn't exist in the array of V0s or was used!" << endl;
2212
2213 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2214 cascade->GetPosCov(covVtx);
2215
2216 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2217 covVtx,
2218 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2219 vcascade,
2220 AliAODVertex::kV0);
2221 vcascade->AddDaughter(vV0FromCascade);
2222 }
2223
2224 // Add the positive tracks from the V0
2225
2226 if (! usedTrack[posFromV0]) {
2227
2228 usedTrack[posFromV0] = kTRUE;
2229
2230 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2231 esdTrack->GetPxPyPz(p);
2232 esdTrack->GetXYZ(pos);
2233 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2234 esdTrack->GetESDpid(pid);
2235
2236 vV0FromCascade->AddDaughter(aodTrack =
2237 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2238 esdTrack->GetLabel(),
2239 p,
2240 kTRUE,
2241 pos,
2242 kFALSE,
2243 covTr,
2244 (Short_t)esdTrack->GetSign(),
2245 esdTrack->GetITSClusterMap(),
2246 pid,
2247 vV0FromCascade,
2248 kTRUE, // check if this is right
2249 kFALSE, // check if this is right
2250 AliAODTrack::kSecondary)
2251 );
2252 aodTrack->ConvertAliPIDtoAODPID();
2253 }
2254 else {
2255 cerr << "Error: event " << iEvent << " cascade " << nCascade
2256 << " track " << posFromV0 << " has already been used!" << endl;
2257 }
2258
2259 // Add the negative tracks from the V0
2260
2261 if (!usedTrack[negFromV0]) {
2262
2263 usedTrack[negFromV0] = kTRUE;
2264
2265 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2266 esdTrack->GetPxPyPz(p);
2267 esdTrack->GetXYZ(pos);
2268 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2269 esdTrack->GetESDpid(pid);
2270
2271 vV0FromCascade->AddDaughter(aodTrack =
2272 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2273 esdTrack->GetLabel(),
2274 p,
2275 kTRUE,
2276 pos,
2277 kFALSE,
2278 covTr,
2279 (Short_t)esdTrack->GetSign(),
2280 esdTrack->GetITSClusterMap(),
2281 pid,
2282 vV0FromCascade,
2283 kTRUE, // check if this is right
2284 kFALSE, // check if this is right
2285 AliAODTrack::kSecondary)
2286 );
2287 aodTrack->ConvertAliPIDtoAODPID();
2288 }
2289 else {
2290 cerr << "Error: event " << iEvent << " cascade " << nCascade
2291 << " track " << negFromV0 << " has already been used!" << endl;
2292 }
2293
2294 // Add the bachelor track from the cascade
2295
2296 Int_t bachelor = cascade->GetBindex();
2297
2298 if(!usedTrack[bachelor]) {
2299
2300 usedTrack[bachelor] = kTRUE;
2301
2302 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2303 esdTrack->GetPxPyPz(p);
2304 esdTrack->GetXYZ(pos);
2305 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2306 esdTrack->GetESDpid(pid);
2307
2308 vcascade->AddDaughter(aodTrack =
2309 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2310 esdTrack->GetLabel(),
2311 p,
2312 kTRUE,
2313 pos,
2314 kFALSE,
2315 covTr,
2316 (Short_t)esdTrack->GetSign(),
2317 esdTrack->GetITSClusterMap(),
2318 pid,
2319 vcascade,
2320 kTRUE, // check if this is right
2321 kFALSE, // check if this is right
2322 AliAODTrack::kSecondary)
2323 );
2324 aodTrack->ConvertAliPIDtoAODPID();
2325 }
2326 else {
2327 cerr << "Error: event " << iEvent << " cascade " << nCascade
2328 << " track " << bachelor << " has already been used!" << endl;
2329 }
2330
2331 // Add the primary track of the cascade (if any)
2332
2333 } // end of the loop on cascades
2334
2335 // V0s
2336
2337 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2338
2339 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2340
2341 AliESDv0 *v0 = esd->GetV0(nV0);
2342
2343 v0->GetXYZ(pos[0], pos[1], pos[2]);
2344 v0->GetPosCov(covVtx);
2345
2346 AliAODVertex * vV0 =
2347 new(vertices[jVertices++]) AliAODVertex(pos,
2348 covVtx,
2349 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2350 primary,
2351 AliAODVertex::kV0);
2352 primary->AddDaughter(vV0);
2353
2354 Int_t posFromV0 = v0->GetPindex();
2355 Int_t negFromV0 = v0->GetNindex();
2356
2357 // Add the positive tracks from the V0
2358
2359 if (!usedTrack[posFromV0]) {
2360
2361 usedTrack[posFromV0] = kTRUE;
2362
2363 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2364 esdTrack->GetPxPyPz(p);
2365 esdTrack->GetXYZ(pos);
2366 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2367 esdTrack->GetESDpid(pid);
2368
2369 vV0->AddDaughter(aodTrack =
2370 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2371 esdTrack->GetLabel(),
2372 p,
2373 kTRUE,
2374 pos,
2375 kFALSE,
2376 covTr,
2377 (Short_t)esdTrack->GetSign(),
2378 esdTrack->GetITSClusterMap(),
2379 pid,
2380 vV0,
2381 kTRUE, // check if this is right
2382 kFALSE, // check if this is right
2383 AliAODTrack::kSecondary)
2384 );
2385 aodTrack->ConvertAliPIDtoAODPID();
2386 }
2387 else {
2388 cerr << "Error: event " << iEvent << " V0 " << nV0
2389 << " track " << posFromV0 << " has already been used!" << endl;
2390 }
2391
2392 // Add the negative tracks from the V0
2393
2394 if (!usedTrack[negFromV0]) {
2395
2396 usedTrack[negFromV0] = kTRUE;
2397
2398 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2399 esdTrack->GetPxPyPz(p);
2400 esdTrack->GetXYZ(pos);
2401 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2402 esdTrack->GetESDpid(pid);
2403
2404 vV0->AddDaughter(aodTrack =
2405 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2406 esdTrack->GetLabel(),
2407 p,
2408 kTRUE,
2409 pos,
2410 kFALSE,
2411 covTr,
2412 (Short_t)esdTrack->GetSign(),
2413 esdTrack->GetITSClusterMap(),
2414 pid,
2415 vV0,
2416 kTRUE, // check if this is right
2417 kFALSE, // check if this is right
2418 AliAODTrack::kSecondary)
2419 );
2420 aodTrack->ConvertAliPIDtoAODPID();
2421 }
2422 else {
2423 cerr << "Error: event " << iEvent << " V0 " << nV0
2424 << " track " << negFromV0 << " has already been used!" << endl;
2425 }
2426
2427 } // end of the loop on V0s
2428
2429 // Kinks: it is a big mess the access to the information in the kinks
2430 // The loop is on the tracks in order to find the mother and daugther of each kink
2431
2432
2433 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2434
2435
2436 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2437
2438 Int_t ikink = esdTrack->GetKinkIndex(0);
2439
2440 if (ikink) {
2441 // Negative kink index: mother, positive: daughter
2442
2443 // Search for the second track of the kink
2444
2445 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2446
2447 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2448
2449 Int_t jkink = esdTrack1->GetKinkIndex(0);
2450
2451 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2452
2453 // The two tracks are from the same kink
2454
2455 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2456
2457 Int_t imother = -1;
2458 Int_t idaughter = -1;
2459
2460 if (ikink<0 && jkink>0) {
2461
2462 imother = iTrack;
2463 idaughter = jTrack;
2464 }
2465 else if (ikink>0 && jkink<0) {
2466
2467 imother = jTrack;
2468 idaughter = iTrack;
2469 }
2470 else {
2471 cerr << "Error: Wrong combination of kink indexes: "
2472 << ikink << " " << jkink << endl;
2473 continue;
2474 }
2475
2476 // Add the mother track
2477
2478 AliAODTrack * mother = NULL;
2479
2480 if (!usedTrack[imother]) {
2481
2482 usedTrack[imother] = kTRUE;
2483
2484 AliESDtrack *esdTrack = esd->GetTrack(imother);
2485 esdTrack->GetPxPyPz(p);
2486 esdTrack->GetXYZ(pos);
2487 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2488 esdTrack->GetESDpid(pid);
2489
2490 mother =
2491 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2492 esdTrack->GetLabel(),
2493 p,
2494 kTRUE,
2495 pos,
2496 kFALSE,
2497 covTr,
2498 (Short_t)esdTrack->GetSign(),
2499 esdTrack->GetITSClusterMap(),
2500 pid,
2501 primary,
2502 kTRUE, // check if this is right
2503 kTRUE, // check if this is right
2504 AliAODTrack::kPrimary);
2505 primary->AddDaughter(mother);
2506 mother->ConvertAliPIDtoAODPID();
2507 }
2508 else {
2509 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2510 << " track " << imother << " has already been used!" << endl;
2511 }
2512
2513 // Add the kink vertex
2514 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2515
2516 AliAODVertex * vkink =
2517 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2518 NULL,
2519 0.,
2520 mother,
2521 AliAODVertex::kKink);
2522 // Add the daughter track
2523
2524 AliAODTrack * daughter = NULL;
2525
2526 if (!usedTrack[idaughter]) {
2527
2528 usedTrack[idaughter] = kTRUE;
2529
2530 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2531 esdTrack->GetPxPyPz(p);
2532 esdTrack->GetXYZ(pos);
2533 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2534 esdTrack->GetESDpid(pid);
2535
2536 daughter =
2537 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2538 esdTrack->GetLabel(),
2539 p,
2540 kTRUE,
2541 pos,
2542 kFALSE,
2543 covTr,
2544 (Short_t)esdTrack->GetSign(),
2545 esdTrack->GetITSClusterMap(),
2546 pid,
2547 vkink,
2548 kTRUE, // check if this is right
2549 kTRUE, // check if this is right
2550 AliAODTrack::kPrimary);
2551 vkink->AddDaughter(daughter);
2552 daughter->ConvertAliPIDtoAODPID();
2553 }
2554 else {
2555 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2556 << " track " << idaughter << " has already been used!" << endl;
2557 }
2558
2559
2560 }
2561 }
2562
2563 }
2564
2565 }
2566
2567
2568 // Tracks (primary and orphan)
2569
2570 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2571
2572
2573 if (usedTrack[nTrack]) continue;
2574
2575 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2576 esdTrack->GetPxPyPz(p);
2577 esdTrack->GetXYZ(pos);
2578 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2579 esdTrack->GetESDpid(pid);
2580
2581 Float_t impactXY, impactZ;
2582
2583 esdTrack->GetImpactParameters(impactXY,impactZ);
2584
2585 if (impactXY<3) {
2586 // track inside the beam pipe
2587
2588 primary->AddDaughter(aodTrack =
2589 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2590 esdTrack->GetLabel(),
2591 p,
2592 kTRUE,
2593 pos,
2594 kFALSE,
2595 covTr,
2596 (Short_t)esdTrack->GetSign(),
2597 esdTrack->GetITSClusterMap(),
2598 pid,
2599 primary,
2600 kTRUE, // check if this is right
2601 kTRUE, // check if this is right
2602 AliAODTrack::kPrimary)
2603 );
2604 aodTrack->ConvertAliPIDtoAODPID();
2605 }
2606 else {
2607 // outside the beam pipe: orphan track
2608 aodTrack =
2609 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2610 esdTrack->GetLabel(),
2611 p,
2612 kTRUE,
2613 pos,
2614 kFALSE,
2615 covTr,
2616 (Short_t)esdTrack->GetSign(),
2617 esdTrack->GetITSClusterMap(),
2618 pid,
2619 NULL,
2620 kFALSE, // check if this is right
2621 kFALSE, // check if this is right
2622 AliAODTrack::kOrphan);
2623 aodTrack->ConvertAliPIDtoAODPID();
2624 }
2625 } // end of loop on tracks
2626
2627 // muon tracks
2628 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2629 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2630
2631 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2632 p[0] = esdMuTrack->Px();
2633 p[1] = esdMuTrack->Py();
2634 p[2] = esdMuTrack->Pz();
2635 pos[0] = primary->GetX();
2636 pos[1] = primary->GetY();
2637 pos[2] = primary->GetZ();
2638
2639 // has to be changed once the muon pid is provided by the ESD
2640 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2641
2642 primary->AddDaughter(
2643 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2644 0, // no label provided
2645 p,
2646 kTRUE,
2647 pos,
2648 kFALSE,
2649 NULL, // no covariance matrix provided
2650 (Short_t)-99, // no charge provided
2651 0, // no ITSClusterMap
2652 pid,
2653 primary,
2654 kTRUE, // check if this is right
2655 kTRUE, // not used for vertex fit
2656 AliAODTrack::kPrimary)
2657 );
2658 }
2659
2660 // Access to the AOD container of clusters
2661 TClonesArray &clusters = *(aod->GetClusters());
2662 Int_t jClusters=0;
2663
2664 // Calo Clusters
2665 Int_t nClusters = esd->GetNumberOfCaloClusters();
2666
2667 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2668
2669 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2670
2671 Int_t id = cluster->GetID();
2672 Int_t label = -1;
2673 Float_t energy = cluster->E();
2674 cluster->GetPosition(posF);
2675 AliAODVertex *prodVertex = primary;
2676 AliAODTrack *primTrack = NULL;
2677 Char_t ttype=AliAODCluster::kUndef;
2678
2679 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2680 else if (cluster->IsEMCAL()) {
2681
2682 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2683 ttype = AliAODCluster::kEMCALPseudoCluster;
2684 else
2685 ttype = AliAODCluster::kEMCALClusterv1;
2686
2687 }
2688
2689 new(clusters[jClusters++]) AliAODCluster(id,
2690 label,
2691 energy,
2692 pos,
2693 NULL, // no covariance matrix provided
2694 NULL, // no pid for clusters provided
2695 prodVertex,
2696 primTrack,
2697 ttype);
2698
2699 } // end of loop on calo clusters
2700
2701 delete [] usedTrack;
2702 delete [] usedV0;
2703 delete [] usedKink;
2704
2705 // fill the tree for this event
2706 aodTree->Fill();
2707 } // end of event loop
2708
2709 aodTree->GetUserInfo()->Add(aod);
2710
2711 // close ESD file
2712 esdFile->Close();
2713
2714 // write the tree to the specified file
2715 aodFile = aodTree->GetCurrentFile();
2716 aodFile->cd();
2717 aodTree->Write();
2718
2719 return;
2720}
2721
2722void AliReconstruction::WriteAlignmentData(AliESD* esd)
2723{
2724 // Write space-points which are then used in the alignment procedures
2725 // For the moment only ITS, TRD and TPC
2726
2727 // Load TOF clusters
2728 if (fTracker[3]){
2729 fLoader[3]->LoadRecPoints("read");
2730 TTree* tree = fLoader[3]->TreeR();
2731 if (!tree) {
2732 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2733 return;
2734 }
2735 fTracker[3]->LoadClusters(tree);
2736 }
2737 Int_t ntracks = esd->GetNumberOfTracks();
2738 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2739 {
2740 AliESDtrack *track = esd->GetTrack(itrack);
2741 Int_t nsp = 0;
2742 Int_t idx[200];
2743 for (Int_t iDet = 3; iDet >= 0; iDet--)
2744 nsp += track->GetNcls(iDet);
2745 if (nsp) {
2746 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2747 track->SetTrackPointArray(sp);
2748 Int_t isptrack = 0;
2749 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2750 AliTracker *tracker = fTracker[iDet];
2751 if (!tracker) continue;
2752 Int_t nspdet = track->GetNcls(iDet);
2753 if (nspdet <= 0) continue;
2754 track->GetClusters(iDet,idx);
2755 AliTrackPoint p;
2756 Int_t isp = 0;
2757 Int_t isp2 = 0;
2758 while (isp < nspdet) {
2759 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2760 const Int_t kNTPCmax = 159;
2761 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
2762 if (!isvalid) continue;
2763 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2764 }
2765 }
2766 }
2767 }
2768 if (fTracker[3]){
2769 fTracker[3]->UnloadClusters();
2770 fLoader[3]->UnloadRecPoints();
2771 }
2772}
2773
2774//_____________________________________________________________________________
2775void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2776{
2777 // The method reads the raw-data error log
2778 // accumulated within the rawReader.
2779 // It extracts the raw-data errors related to
2780 // the current event and stores them into
2781 // a TClonesArray inside the esd object.
2782
2783 if (!fRawReader) return;
2784
2785 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2786
2787 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2788 if (!log) continue;
2789 if (iEvent != log->GetEventNumber()) continue;
2790
2791 esd->AddRawDataErrorLog(log);
2792 }
2793
2794}
2795
2796TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2797 // Dump a file content into a char in TNamed
2798 ifstream in;
2799 in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2800 Int_t kBytes = (Int_t)in.tellg();
2801 printf("Size: %d \n",kBytes);
2802 TNamed *fn = 0;
2803 if(in.good()){
2804 char* memblock = new char [kBytes];
2805 in.seekg (0, ios::beg);
2806 in.read (memblock, kBytes);
2807 in.close();
2808 TString fData(memblock,kBytes);
2809 fn = new TNamed(fName,fData);
2810 printf("fData Size: %d \n",fData.Sizeof());
2811 printf("fName Size: %d \n",fName.Sizeof());
2812 printf("fn Size: %d \n",fn->Sizeof());
2813 delete[] memblock;
2814 }
2815 else{
2816 AliInfo(Form("Could not Open %s\n",fPath.Data()));
2817 }
2818
2819 return fn;
2820}
2821
2822void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2823 // This is not really needed in AliReconstruction at the moment
2824 // but can serve as a template
2825
2826 TList *fList = fTree->GetUserInfo();
2827 TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2828 printf("fn Size: %d \n",fn->Sizeof());
2829
2830 TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2831 const char* cdata = fn->GetTitle();
2832 printf("fTmp Size %d\n",fTmp.Sizeof());
2833
2834 int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2835 printf("calculated size %d\n",size);
2836 ofstream out(fName.Data(),ios::out | ios::binary);
2837 out.write(cdata,size);
2838 out.close();
2839
2840}
2841