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